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Abstract 

A  methodology  is  developed  for  determining  the  validity  of  making  a  statistical 
turbulent  approach  using  Kolmogorov  theory  to  an  aero-optical  turbulent  flow.  Kol¬ 
mogorov  theory  provides  a  stochastic  method  that  has  a  greatly  simplified  and  robust 
method  for  calculating  atmospheric  turbulence  effects  on  optical  beam  propagation, 
which  could  simplify  similar  approaches  to  chaotic  aero-optical  flows.  A  2-D  lami¬ 
nar  Navier-Stokes  CFD  Solver  (AVUS)  is  run  over  a  splitter  plate  type  geometry  to 
create  an  aero-optical  like  shear  mixing  layer  turbulence  held.  A  Matlab  algorithm 
is  developed  to  import  the  how  data  and  calculates  the  structure  functions,  struc¬ 
ture  constant,  and  Fried  Parameter  (r0)  and  compares  them  to  expected  Kolmogorov 
distributions  assuming  an  r2//3  power  law.  The  range  of  C2’ s  developed  from  the 
structure  functions  are  not  constant  with  separation  distance,  and  ranged  between 
l(T12-l(r10.  There  is  a  consistent  range  of  data  overlap  within  the  C2’ s  derived  from 
various  methods  for  separation  distances  within  the  range  0.01m-0.02m.  Within  this 
range  rQ  is  found  to  be  approximately  0.05m  which  is  a  reasonable  value.  This  par¬ 
ticular  2-D  shear  mixing  layer  was  found  to  be  non-Kolmogorov,  but  further  grid 
refinement  and  data  sampling  may  provide  a  more  Kolmogorov  like  distribution. 
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NUMERICAL  INVESTIGATION  OF  STATISTICAL  TURBULENCE  EFFECTS 
ON  BEAM  PROPAGATION  THROUGH  2-D  SHEAR  MIXING  LAYER 

I.  Introduction 


1.1  Research  Overview 

The  advent  of  laser  technology  was  described  as  a  solution  without  a  problem. 
The  Air  Force  has  made  significant  investments  into  research  of  new  applications 
for  laser  technologies  in  meeting  its  day  to  day  mission  requirements.  Researchers 
are  investigating  ways  to  use  lasers  for  air  and  space  communication,  asset  defense, 
and  airborne  high  power  directed  energy  effects.  The  critical  aspect  in  using  lasers 
in  this  manner  is  that  they  involve  propagating  the  beam  through  the  atmosphere, 
which  is  not  necessarily  the  most  accommodating  medium.  The  atmosphere  has  many 
deleterious  effects  on  laser  beam  propagation.  These  deletrious  effects  arise  from  the 
spectral  wavelength  of  the  source  and  changes  in  the  index  of  refraction  from  density 
changes  in  the  medium.  They  are  directly  attributed  to  three  physical  phenomenon: 
absorption,  scattering,  and  atmospheric  turbulence.  These  effects  prevent  the  laser 
from  putting  the  entirety  of  its  energy  or  radiant  intensity  on  a  target  or  cause 
the  beam  to  be  directed  away  from  the  target,  which  greatly  reduces  the  ability  to 
precisely  distribute  and  introduce  the  effects  desired.  By  understanding  how  these 
mechanisms  affect  propagation,  these  beam  quality  issues  can  be  better  anticipated 
and  mitigated. 

While  absorption  and  scattering  are  important  phenomenon,  the  focus  of  this 
research  is  turbulence.  Turbulence  is  instability  in  motion  of  fluid,  or  the  atmosphere. 
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Turbulence  arises  from  injection  of  energy  into  the  fluid  causing  the  motion  to  become 
unstable.  This  source  of  this  energy  injection  is  usually  introduced  thermally  from 
temperature  gradients  with  height  in  the  atmosphere  or  aerodynamic  friction  along 
the  ground,  ocean,  the  skin  of  an  aircraft,  or  any  other  surface.  Turbulence  is  often 
represented  physically  by  whorl-like  structures  called  eddies.  These  eddies  come  in 
many  shapes  and  sizes  which  leads  to  a  need  to  describe  these  eddies  with  respect  to 
a  scale.  In  the  instances  above,  atmospheric  eddies  within  the  inertial  sub  range  can 
range  from  200  meters  to  0.002  meters  compared  to  the  aerodynamic  eddies,  which 
can  be  described  as  microscopic  as  they  are  on  the  order  of  2  millimeters  and  smaller. 

Being  able  to  model  and  predict  turbulence  is  quite  challenging  because  of  its 
large  range  of  length  and  time  scales.  The  injection  of  kinetic  energy  to  produce  the 
larger  eddies  begins  to  fade  into  to  smaller  and  smaller  eddies.  Because  of  the  small 
size  scale  of  the  initial  eddies  and  the  time  scale  over  which  the  energy  is  injected  into 
the  fluid  makes  it  difficult  to  describe  the  behavior  numerically  real-time.  Computing 
turbulence  real-time  would  require  a  very  fine  resolution  evolution  in  time  and  space, 
and  would  only  be  accurate  for  a  very  small  period  of  time.  This  would  require 
a  large  number  of  computational  resources  and  time  to  compute.  There  are  other 
numerical  methods,  such  as  Large  Eddy  Simulation,  that  model  turbulence  behavior 
as  chaotic  or  random  while  the  fluid  flow  evolves  as  solved  using  the  Navier-Stokes 
equations.  Therefore,  statistical  or  stochastic  methods  such  as  Kolmogorov’s  method 
were  derived  for  the  purpose  of  simulating  beam  propagation  effects. 

The  heart  of  this  research  is  investigating  the  effects  the  resultant  turbulence  fields 
from  these  different  modeling  approaches  has  on  propagating  a  laser  through  the  at¬ 
mosphere.  Does  micro-scale  aerodynamic  turbulence  described  chaotically  produce 
turbulence  fields  that  affect  simulated  laser  beam  propagation  in  a  quantifiably  dif¬ 
ferent  way  than  atmospheric  turbulence  modeled  statistically?  By  classifying  which 
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model  has  the  more  dominant  effect  on  the  beam  could  change  the  approach  engi¬ 
neers  take  to  mitigate  atmospheric  turbulence  effects  on  propagation.  It  may  also 
cause  researchers  to  reevaluate  approaches  to  analyzing  the  modeling  descriptions  of 
atmospheric  turbulence. 

1.2  Brief  Background 

1.2.1  Atmospheric  Effects  on  Lasers. 

Laser  is  the  short  hand  form  for  Light  Amplification  through  Stimulated  Emis¬ 
sion  of  Radiation.  A  laser  works  by  taking  an  input  electromagnetic  radiation  source 
and  amplifies  it  by  producing  additional  electromagnetic  radiation  produced  by  the 
quantum  mechanical  means  of  stimulated  emission  through  a  receptive  gain  medium. 
The  major  aspect  of  this  process  is  that  it  leads  to  an  increased  irradiance  of  electro¬ 
magnetic  radiation  which  carries  increased  spectral  power  in  a  specific  wavelength, 
frequency,  and  polarization.  Notice  the  word  light  has  not  been  used  in  the  above  de¬ 
scription.  This  is  to  keep  the  discussion  relevant  to  any  radiation  that  falls  within  the 
electromagnetic  spectrum.  The  nature  of  this  electromagnetic  radiation  in  regards  to 
light  is  rather  complex  with  regard  to  it  being  made  up  of  corpuscular  components 
referred  to  as  photons  or  as  waves.  Regardless  of  which  form  actually  describes  light  it 
is  important  to  understand  that  this  amplified  radiation  is  susceptible  to  phenomenon 
that  can  alter  the  amount  of  radiation  that  arrives  to  its  final  destination  or  where 
that  final  destination  is. 

Diffraction  is  the  redirection  of  the  wave  front  or  ray  due  to  a  change  in  the  optical 
properties  of  the  medium  through  which  the  wave  or  ray  is  traveling.  The  optical 
property  of  the  medium  is  often  referred  to  as  the  index  of  refraction.  The  refractive 
index  is  typically  dependent  on  wavelength  of  the  incident  radiation  and  the  density 
of  the  medium.  Therefore,  slight  density  changes  in,  for  example,  the  atmosphere 
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due  to  turbulence,  affect  the  direction  and  intensity  of  the  beam  as  it  traverses  the 
medium. 

Atmospheric  extinction  is  the  reduction  of  the  incident  radiation  through  the 
atmosphere  due  to  absorption  and  scattering  effects.  In  the  same  manner  that  the 
laser  beam  is  created  by  the  quantum  mechanical  means  of  stimulated  emission,  it 
is  also  susceptible  to  giving  up  its  energy  to  molecules  along  the  path  that  would 
be  receptive  to  the  quantum  mechanism  of  absorption.  The  energy  can  also  be  lost 
through  absorption  by  inducing  motion  of  the  molecule  through  rotation  or  vibration. 
This  is  a  major  consideration  due  to  the  number  of  different  molecules  that  constitute 
air,  as  well  as  water  and  the  numerous  available  aerosols  in  the  atmosphere.  The 
nature  of  these  losses  is  primarily  wavelength  driven  and  is  usually  mitigated  in  the 
operational  design  wavelength  bandwidth  of  the  laser.  [11] 

Absorption  also  contributes  to  another  deleterious  effect  on  beam  propagation  in 
the  form  of  thermal  blooming.  As  the  air  and  other  atmospheric  constituents  absorb 
the  incident  radiation,  that  energy  often  takes  the  form  of  heat.  Heat  produces 
thermal  gradients  in  the  air  which  results  in  density  changes  in  the  medium.  As 
stated  above,  these  density  changes  result  in  local  changes  in  the  index  of  refraction, 
which  can  shift  the  direction  of  the  beam.  [11] 

Scattering  is  an  extinction  mechanism  in  that  it  removes  energy  from  the  laser 
beam  by  redirecting  the  radiation  through  elastic  collisions  with  present  particles 
in  directions  not  in  the  direction  of  propagation.  While  the  incident  radiation  is  not 
necessarily  lost  or  destroyed  in  the  action,  the  beam  has  now  lost  some  of  its  constitu¬ 
tional  integrity,  reducing  the  amount  of  radiation  on  target.  This  mechanism  is  again 
primarily  wavelength  dependent  and  can  be  compensated  by  design  considerations 
for  the  operational  window  of  the  laser.  [11] 

These  mechanisms  represent  a  few  of  the  design  challenges  facing  engineers  in 
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the  incorporation  of  laser  systems  for  military  application  in  air  and  space.  The 
above  effects  are  losses  in  energy  due  to  atmospheric  extinction  which  are  mostly 
wavelength  dependent.  However,  diffraction  due  to  refractive  index  changes  is  driven 
by  density  changes.  Turbulence  must  be  carefully  considered  in  system  design,  as 
turbulence  introduces  varying  local  density  changes,  providing  a  dynamic  refractive 
index  change.  By  better  understanding  sources  of  turbulence  and  how  to  describe 
them  are  crucial  when  trying  to  address  laser  integration  into  Air  Force  systems. 

1.2.2  Turbulence. 

As  mentioned  earlier,  turbulence  is  instability  in  fluid  flow,  which  causes  localized 
variances  in  fluid  speed  and  density.  The  source  of  these  instabilities  is  the  intro¬ 
duction  of  additional  kinetic  energy  into  the  flow.  Atmospheric  turbulence  can  be 
characterized  by  three  means:  mechanical,  thermal,  and  inertial.  Mechanical  turbu¬ 
lence  is  the  result  of  variance  in  the  wind  speed  at  various  layers  referred  to  as  shear. 
The  changes  in  shear  can  be  attributed  to  surface  drag,  wake  flow,  or  free  stream 
shear.  Surface  drag  creates  variations  due  to  slower  wind  speeds  at  the  surface,  due 
to  friction,  than  speeds  higher  aloft.  Wake  flow  is  effect  of  surface  roughness  or 
protuberances  from  the  surfaces,  such  as  trees  or  buildings,  which  cause  wind  speed 
variations  as  the  fluid  must  move  around  these  obstacles.  Free  stream  shear  is  wind 
speed  variations  that  occur  naturally  away  from  any  solid  surface.  Thermal  turbu¬ 
lence  is  result  of  changes  in  temperature,  usually  from  radiative  heating  or  cooling, 
that  cause  air  to  rise  or  sink  buoyantly.  Inertial  turbulence  is  the  product  of  the  loss 
of  kinetic  energy  of  eddies,  which  leads  to  a  dissipation  of  eddies  into  smaller  and 
smaller  eddies.  [20] 

In  light  of  these  mechanisms  for  creating  or  affecting  existing  turbulent  sources, 
turbulence  is  highly  dynamic  in  nature,  and  hard  to  predict.  In  this  regard,  the 
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constantly  changing  size  and  time  scales  for  how  the  dissipation  occurs  is  hard  to 
model  and  predict  as  these  changes  are  often  random  or  non-linear,  and  require  very 
small  scale  spatial  resolution  to  capture  the  physics  or  become  very  inaccurate  for 
longer  time  scales.  This  inability  to  accurately  capture  turbulent  flow  behavior  led  to 
two  very  different  approaches  in  modeling  turbulence:  chaotic,  allowing  the  random 
behavior  to  evolve  in  the  flow,  or  statistical,  developing  a  parameterized  description 
of  the  turbulence  with  regard  to  temperature  range  and  average  velocities.  [20] 

The  current  problem  in  approaching  turbulence  effects  on  laser  beam  propagation 
is  that  there  are  two  ranges  or  regions  of  turbulent  flows  that  have  contrasting  view¬ 
points,  scales,  and  assumptions  made  that  has  divided  the  research  community  in  this 
field.  The  first  viewpoint  is  that  atmospheric  turbulence  is  the  dominant  factor.  The 
strength  of  this  approach  is  that  because  of  the  ranges  and  somewhat  predictable 
nature  of  the  atmosphere  allows  certain  assumptions  to  be  made  which  allows  for 
a  more  simplified  or  standardized  approach  to  addressing  turbulent  effects  on  beam 
propagation.  The  weakness  is  that  it  ignores  small  scale  turbulent  features  that  may 
be  present  in,  for  example,  an  airborne  system  that  has  boundary  layer  turbulence 
over  an  aperture.  The  second  viewpoint  is  that  these  aero-optical  features  cannot  be 
ignored.  The  strength  of  this  approach  is  that  it  allows  for  a  more  complete  account 
for  turbulent  effects  on  beam  propagation.  The  problem  though  is  that  the  smaller 
scale  and  unpredictable  nature  of  these  aero-optical  eddies  make  it  difficult  to  make 
the  same  assumptions  that  a  statistical  viewpoint  uses. 

The  presiding  theory  for  the  statistical  representation  of  the  atmosphere  is  the 
Kolmogorov  Theory  of  1941.  Kolmogorov  theory  provides  a  way  to  parameterize  the 
eddy  size  scale,  velocity  and  temperature  profiles,  to  provide  structure  constants  that 
describe  the  nature  of  the  turbulence  within  that  range.  Using  these  structure  terms, 
optical  parameters  can  be  developed  and  used  to  predict  the  behavior  of  a  wave  front 
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to  propagate  through  the  turbulence.  This  theory  is  very  useful;  however  it  is  only 
applicable  to  modeling  turbulence  effects  on  optical  performance  within  the  inertial 
sub-range  of  turbulence.  [2] 

A  lot  of  the  analytical  work  for  investigating  aero-optical  turbulence  is  performed 
using  numerical  models.  Modeling  numerical  fluid  flow  can  be  accomplished  using 
the  Navier-Stokes  equations.  Numerical  modeling  of  fluid  dynamics  using  the  Navier- 
Stokes  equations  has  become  very  robust,  and  there  are  many  commercially  available 
flow  solvers  that  allow  these  flow  physics  to  be  calculated.  As  understanding  of  tur¬ 
bulence  behavior  has  evolved,  mathematical  models  that  predict  the  unsteady  nature 
of  turbulence  have  been  incorporated  with  these  Navier-Stokes  solvers  to  provide  a 
characteristic  representation  of  the  random  nature  of  turbulence.  These  computa¬ 
tional  flow  simulations  provide  some  opportunities  to  investigate  the  fluid  dynamics 
around  an  airborne  system.  As  previously  stated  though,  these  simulations  can  be 
very  time  consuming  to  generate  as  they  require  small  temporal  and  spatial  scales, 
and  large  amounts  of  processing  to  generate.  Even  then  these  simulations  may  only 
be  qualitatively  representative  rather  than  quantitatively  useful.  Therefore,  these 
simulations  are  representations  of  representative  flow  condition  or  used  in  diagnostic 
of  what  the  flow  may  be  like,  and  not  neccesarily  useful  for  a  rigorous  or  accurate 
accessment  of  turbulent  effects. 

However,  if  there  is  little  difference  between  a  characteristic  aero-optical  flowfield 
and  its  statistical  representation  then  there  would  be  an  opportunity  to  reconcile 
these  two  approaches  and  simplify  the  way  engineers  design  and  compensate  for  both 
atmospheric  and  aero-optical  turbulence.  This  is  fundamental  motivation  for  this 
research. 
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1.3  Investigative  Setup  and  Scope 


In  order  to  create  the  basis  for  this  comparison  of  turbulence  effects,  a  fluid 
flow  similar  to  an  aero-optic  flow  will  be  generated  using  Laminar  Navier-Stokes 
approaches  that  are  available  in  Computational  Fluid  Dynamics  (CFD)  software.  A 
two-dimensional  (2-D)  splitter  plate  type  geometry  to  generate  a  laminar  shear  layer 
as  representative  turbulent  field  for  analysis.  The  output  from  the  fluid  model  will 
then  be  accessed  using  stochastic  methodologies  to  determine  if  this  characteristic 
aero-optical  turbulent  flow  field  has  Kolmogorov  like  qualities. 

1.4  Outline  of  Thesis 

The  remainder  of  this  document  is  divided  among  4  chapters  following  this  intro¬ 
duction.  Chapter  two  is  an  in  depth  discussion  of  all  relevant  theory  to  the  investi¬ 
gation,  as  well  as  a  look  at  previous  or  current  research  into  this  particular  approach 
or  topic.  Chapter  three  will  state  with  more  detail  the  research  approach,  provid¬ 
ing  a  more  detailed  technical  description  of  software  used,  and  specific  parameters 
being  considered  for  matters  of  comparison.  Chapter  four  will  present  the  results 
of  the  conducted  research  and  will  provide  a  thorough  analysis  of  the  collected  data 
with  regards  to  the  nature  of  the  research.  Chapter  five  will  conclude  this  thesis.  It 
will  provide  a  short  summary  of  results,  and  provide  insights  with  regards  to  future 
research  along  this  topic. 


II.  Background  and  Theory 


This  chapter  presents  the  reader  with  a  review  of  any  pertinent  background,  the¬ 
ory,  and  research  previously  conducted  in  this  field  of  study.  Chapter  two  begins 
with  a  discussion  of  atmospheric  effects  on  laser  beam  propagation.  A  more  in  depth 
discussion  of  turbulence  is  then  provided.  After  establishing  a  background  on  turbu¬ 
lence,  a  review  of  Kolmogorov’s  Theory  is  provided  as  means  to  discuss  a  statistical 
means  for  analyzing  turbulence  effects  on  beam  propagation.  A  summary  is  then  pro¬ 
vided  presenting  researched  examples  of  where  Kolmogorov’s  theory  is  not  capable  of 
capturing  the  effects  of  turbulence  on  beam  propagation.  There  is  then  a  discussion 
of  fluid  mechanical  studies  on  turbulence  effects  for  aero-optical  applications.  The 
chapter  finishes  with  a  quick  discussion  of  Computational  Fluid  Dynamics  (CFD) 
simulation  tools  to  be  used  in  this  study. 

2.1  Atmospheric  Effects  on  Lasers 

Light  is  susceptible  to  many  forms  of  alteration  due  to  diffraction  and  attenuation. 
The  atmosphere  can  produce  the  means  to  introduce  these  influences  in  the  form  of 
atmospheric  extinction  due  to  scatterring  and  absorption,  and  atmospheric  turbu¬ 
lence.  These  mechanisms  represent  a  few  of  the  design  challenges  facing  engineers  in 
the  incorporation  of  laser  systems  for  military  application  in  air  and  space.  Atmo¬ 
spheric  extinction  effects  are  mostly  wavelength  dependent.  Atmospheric  turbulence 
introduces  density  changes  in  the  medium  leading  to  diffraction  from  refractive  index 
changes.  An  example  of  how  optical  density  affects  wave  propagation  can  be  seen  in 
Figures  1  and  2.  As  optical  density  increases  the  transmitted  wavefronts  rotate  away 
from  the  incident  axis  as  seen  in  Figure  1.  With  the  dynamic  nature  of  turbulence, 
the  optical  density  gradient  varies  with  time,  which  creates  transmission  fluctuations 
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as  seen  in  Figure  2.  Turbulence  influence  on  beam  performance  is  difficult  to  mitigate 
because  of  its  inherently  unpredictable  nature.  By  better  understanding  sources  of 
turbulence  and  how  to  describe  them  are  crucial  when  trying  to  address  laser  system 
integration  for  Air  Force  applications.  [11] 

Incident  Waves  Less  Optically 

Dense  Region  Transmitted  Waves 


Figure  1.  Illustration  of  Non-Dynamic  Refractive  Medium  on  Propagation 


Turbulent  Induced 

Incident  Waves  Optical  Density  Changes  Transmitted  Waves 


Figure  2.  Illustration  of  Time  Varying  Refractive  Medium  on  Propagation 
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2.2  Turbulence 


In  order  to  analyze  the  effects  of  turbulence  on  beam  propagation,  it  is  important 
to  understand  what  turbulence  is,  where  it  comes  from,  and  some  representative 
characteristics  that  represent  it.  Turbulence  derives  from  instabilities  in  fluid  flow. 
It  is  characterized  by  irregularity,  or  non-linear  behavior.  This  randomness  results 
in  the  development  of  rotational  structure  called  eddies,  swirling  in  all  directions. 
These  whirling  eddies  make  the  turbulent  flow  diffusive,  causing  rapid  mixing  which 
increases  the  rates  of  transfer  in  momentum,  heat,  and  mass.  Turbulence  is  also 
dissipative  in  that  viscous  forces  increase  the  energy  of  the  fluid,  which  causes  kinetic 
energy  of  the  turbulence  to  decrease  in  order  to  maintain  energy  conservation.  If 
there  is  no  further  injection  of  kinetic  energy  into  the  flow,  then  turbulence  will  decay 
rapidly.  [17] 

It  is  now  important  to  understand  what  causes  turbulence.  The  source  of  at¬ 
mospheric  turbulence  is  the  introduction  of  additional  kinetic  energy  into  the  flow. 
Turbulence  can  be  introduced  by  three  means:  mechanical,  thermal,  and  inertial. 
Mechanical  turbulence  is  the  result  of  variance  in  the  flow  speed  at  various  layers 
referred  to  as  shear.  The  changes  in  shear  can  be  attributed  to  surface  drag,  wake 
flow,  or  free  stream  shear.  In  atmospheric  turbulence,  surface  drag  creates  variations 
due  to  slower  wind  speeds  at  the  surface,  due  to  friction,  than  speeds  higher  aloft. 
Wake  flow  is  effect  of  surface  roughness  or  protuberances  from  the  surfaces,  which 
cause  wind  speed  variations  as  the  fluid  must  move  around  these  obstacles.  Free 
stream  shear  is  just  level  wind  speed  variations  that  occur  naturally  away  from  any 
solid  surface.  Thermal  turbulence  is  result  of  changes  in  temperature,  usually  from 
radiative  heating  or  cooling  in  the  atmosphere,  that  cause  air  to  rise  or  sink  buoy¬ 
antly.  Inertial  turbulence  is  the  product  of  the  loss  of  kinetic  energy  of  eddies,  which 
leads  to  a  dissipation  of  eddies  into  smaller  and  smaller  eddies  as  describe  above.  [20] 
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These  characteristic  make  it  easy  to  see  that  turbulence  by  nature  is  very  non¬ 
linear.  This  makes  it  very  difficult  if  not  impossible  to  describe  turbulence  spatially  or 
temporally  with  any  precision.  This  inability  to  accurately  capture  turbulent  flow  be¬ 
havior  led  to  two  very  different  approaches  in  modeling  turbulence:  a  chaotic  aproach 
which  allows  the  random  flow  behavior  to  evolve  using  some  numerical  microphysics 
model,  or  statistically  developing  a  parameterized  description  of  the  turbulence  with 
regard  to  temperature  range,  eddy  scales,  and  average  velocities. 

Turbulence  was  first  systematically  investigated  by  Osborne  Reynolds  in  1883. 
His  major  contribution  to  the  field  was  the  introduction  of  a  nondimensional  means 
for  characterizing  turbulence. 

ul  Inertial  Transport 

Re  —  —  =  — - — -  (1) 

v  Viscous  Transport 

The  Reynolds  Number  (Re)  as  seen  in  Equation  1  is  determined  by  the  flow 

speed  (u),  scale  length  (/),  and  viscosity  (v).  The  Reynolds  number  can  provide 
a  classification  of  the  strength  of  turbulence.  Lower  numbers  represent  a  laminar 
turbulent  flow  which  is  more  smooth  and  less  irregular.  Higher  numbers  represent 
strong  turbulence  characterized  by  random  and  irregular  flow.  [4] 

As  the  Reynolds  number  increases,  the  separation  of  length  and  time  scales  in¬ 
creases  as  well.  These  increases  introduce  more  chaotic  flow  conditions  with  them 
that  are  hard  to  define  mathematically.  In  order  to  begin  a  statistical  assessment 
of  turbulence,  a  dimensional  analysis  is  necessary.  This  involves  the  development  of 
scales  for  the  size  of  the  turbulent  eddies.  The  eddies  which  constitute  turbulence 
represent  a  spectrum  as  seen  in  Figure  3.  This  spectrum  can  be  classified  by  three 
regions:  The  Production  Subrange,  the  Inertial  Subrange,  and  the  Dissipation  Sub¬ 
range.  The  production  subrange  is  where  kinetic  energy  is  introduced  in  the  flow  to 
create  turbulent  eddies.  Turbulence  in  this  region  is  very  anisotropic  and  not  gov- 
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erned  by  any  general  description.  Eventually  that  turbulent  kinetic  energy  starts  to 
be  lost  into  the  internal  energy  of  the  flow  by  molecular  viscosity.  When  those  viscous 
forces  dominate  the  kinetic  energy,  the  flow  enters  the  dissipation  subrange.  This  re¬ 
gion  represents  the  smallest  eddy  structures.  In  between  these  regions  is  the  inertial 
subrange.  This  subrange  is  where  the  energy  lost  to  heat  is  very  small  compared  to 
the  transfer  of  energy  into  the  turbulent  scales.  [20;  8] 


Large  Eddies 


EddyScale 


Small  Eddies 


Figure  3.  Cascade  of  Turbulent  Kinetic  Energy 

Now  that  these  subregions  for  the  spectrum  have  been  introduced,  it  is  important 
to  know  where  sizing  scale  is  derived  from.  Figure  4  is  a  representation  of  how  scaling 
size  originated  from  the  cascade  of  eddies  through  the  turbulence  spectrum.  The 
inertial  subrange  begins  when  there  is  no  further  injection  of  kinetic  energy  into  the 
flow.  The  upper  boundary  (La)  represents  the  size  of  the  largest  eddy  present  as 
energy  injection  ceases.  The  kinetic  energy  continues  to  dissipate  from  this  point, 
however  it  is  inviscid,  energy  transfer  into  turbulent  scales  is  greater  than  energy 
transfer  into  heat.  The  lower  bound  (l0)  is  the  smallest  size  that  eddies  will  take 
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within  the  inertial  subrange.  Below  lQ  the  dissipation  mechanism  is  friction  due  to 
viscosity.  The  lower  boundary  is  usually  the  transition  from  laminar  to  turbulent 
flow.  When  the  eddy  size  is  between  these  bounds  the  turbulence  is  considered  to  be 
in  the  inertial  subrange.  The  utility  is  that  theory  suggests  that  by  establishing  these 
boundary  scales,  the  inertial  subrange  is  homogenous  and  isotropic  which  is  useful  in 
developing  a  statistical  representation  for  the  flow  conditions  that  may  affect  beam 
propagation.  [8] 
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Figure  4. 


Depiction  of  Eddy  Dissipation  and  Scales 


2.3  Statistical  Representation  of  Turbulence 
2.3.1  Defining  Optical  Parameters. 

Before  discussing  Kolmogorov  Theory,  it  is  important  to  understand  some  of  the 
optical  parameters  that  are  used  to  determine  the  statistical  representations  of  a 
turbulent  field  and  how  they  are  derived. 
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2.3. 1.1  Index  of  Refraction. 


The  refractive  index  (n)  of  a  medium  is  the  degree  to  which  a  wavefront  will  be 
redirected  as  it  propagates  as  presented  and  illustrated  previously.  For  free  space  or 
a  vacuum  there  is  no  refraction  and  the  refractive  index  is  set  to  unity.  Stable  air  has 
a  refractive  index  that  is  very  close  to  unity.  However  turbulence  and  atmospheric  air 
constituents  provide  variations  in  density  and  wavelength  dependent  affects  that  can 
create  slight  deviations  from  that  of  stable  air,  which  effects  propagation.  There  are 
several  empirically  derived  methods  for  finding  the  local  index  of  refraction  at  a  point 
in  the  atmosphere  depending  on  beam  wavelength,  and  local  pressure,  temperature, 
and  density.  The  first  is  the  Ecllen  formula  as  seen  in  Equation  2, 


n(x,y,z,t ) 


1  +  (10-s)(8342.12  + 


2406030 
130-  J 


15997  p(x,y,z,t ) 

38.9  —  p0 


(2) 


which  is  dependent  on  the  wavelength  (A  in  pm)  the  spatial  value  of  density  (p(x,  y ,  z,  t  j) 
and  the  reference  or  standard  density  (p0).  The  second  is  the  Gladstone-Dale  Method 
as  seen  Equation  3, 

n(x,  y,z,t)  =  1  +  KGDp(x ,  y,  z,  t )  (3) 


which  uses  a  constant  ( Kqd )  which  has  been  derived  empirically  and  the  spatial  value 
of  density  (p(x,y,  z,t)).  The  Gladstone-Dale  for  a  1  micron  beam  is  approximately 
2.25  x  10“4^.  [21;  12] 

2. 3. 1.2  Phase,  Optical  Path  Length,  and  Optical  Path  Difference. 

Turbulence  distorts  the  phase  front  of  a  plane  wave.  There  are  two  mathematical 
means  that  can  be  used  to  describe  the  amount  of  phase  (0)  alteration  that  occurs 
due  to  turbulent  abberations.  The  first  is  the  Optical  Path  Length  (OPL).  The  OPL 
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through  a  medium  and  be  approximated  using  Equation  4, 


OPL(x,y,z,t)  ~ 


n(x,  y,  z,  t )  ds 


(4) 


where  s  represents  the  propagation  path  length  and  n(x,  y,  z,  t )  represents  the  spatial 
refractive  index  along  that  path.  The  wave’s  phase  (0)  can  be  found  from  the  OPL 
using  Equation  5, 

(j)  =  k0OPL(x,y,z,t)  (5) 

where  kQ  is  the  wavenumber  (k0  =  ).  The  second  means  of  mathematical  describing 

the  wavefront  distortion  is  the  Optical  Path  Distance  (OPD)  which  is  the  difference 
in  the  OPL  relative  to  the  mean  over  an  apperture  as  seen  below  in  Equation  6, 


OPD(x,  y,  z,  t )  =  OPL(x,  y,  z,  t )  —  OPL(s)  (6) 

where  OPL(s)  represents  the  average  of  the  OPL’s  along  a  simulated  plane  perpen¬ 
dicular  to  the  beam  path.  [19] 

2.3.2  Kolmogorov  Theory. 

At  the  forefront  of  this  statistical  representation  of  Turbulence  is  the  Kolmogorov 
Theory  of  1941.  Kolmogorov  theory  provides  a  way  to  parameterize  the  eddy  size 
scale,  velocity  and  temperature  profiles,  to  provide  structure  constants  that  describe 
the  nature  of  the  turbulence  within  that  range.  The  boundary  for  this  characterization 
is  dependent  on  the  scale  size  of  the  turbulent  eddies  as  described  above.  Kolmogorov 
theory  of  1941  provided  two  definitions  from  which  to  make  the  determination  that  a 
turbulence  distribution  could  be  called  homogeneous  and  isotropic.  These  definitions 
define  the  boundaries  of  the  inertial  subrange.  Kolmogorov  also  contributed  two 
hypothesis  on  similarity  which  provides  a  way  to  develop  parameters  for  analysis 
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assuming  a  homogeneous,  isotropic  distribution.  [9] 

The  first  similarity  hypothesis  says  that  at  a  significantly  high  Reynolds  number 
an  isotropic  turbulence  distribution  can  be  defined  by  the  values  of  e  ,  the  average 
energy  dissipation  rate,  and  v,  the  viscosity.  The  average  energy  dissipation  rate  (e)is 
defined  by  the  mean  fluid  velocity  (u)  and  the  outer  scale  length  (L0)  by  the  following 
relationship  in  Equation  7: 


£  ~ 


(7) 


The  second  similarity  hypothesis  states  that  flows  with  sufficiently  high  Reynold’s 
Numbers,  within  the  inertial  subrange  where  the  dominant  turbulence  process  is 
inertial  transfer  in  to  turbulent  scales  rather  than  in  to  viscous  creation  of  heat,  then 
the  distribution  is  dependent  only  on  the  energy  dissipation  rate  (f).  [10] 

A  representative  example  of  this  can  be  found  in  Figure  5.  The  Energy  Spectrum 
should  look  familiar.  It  is  much  like  that  presented  in  Figure  3.  The  major  differences 
between  these  representation  of  scales  as  Figure  5  is  semilog,  and  that  the  indepen¬ 
dent  axis  is  wavenumber  as  opposed  to  eddy  size.  Using  these  similiarity  hypothesis 
Kolmogorov  was  able  to  build  an  expression  for  the  kinetic  energy  spectral  density 
as  presented  in  Equation  8 

E(ka)  =  Kehf  (8) 


The  exponents  in  Equation  8  are  found  using  dimesional  analysis.  The  K  in  Equa¬ 
tion  8  represents  a  constant.  This  constant  is  what  makes  the  curve  fit  as  seen  in 
Figure  5  to  the  spectral  density  function  within  the  inertial  subrange.  It  is  this  ap¬ 
proach  that  provides  a  method  for  the  development  of  structure  functions  and  the 
associated  fit  parameters  such  as  C\.  [14] 
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Figure  5.  Turbulent  Energy  Spectrum  and  Fit  Example.  [14] 

2.3. 2.1  Structure  Functions. 

A  structure  function  is  a  statistical  discription  of  the  average  of  the  differences 
squared  of  two  parameters  separated  by  a  distance,  r.  Equation  9,  which  is  an  auto¬ 
correlation  function  commonly  used  in  turbulence  studies,  provides  a  representative 
view  of  this  concept  where  x  represents  some  parameter  (velocity,  phase,  etc)  and  r 
represents  the  distance  between  the  reference  and  the  point  of  interest. 

Dx(r)  =  (x(0)  -  x(r))2  (9) 

Figure  6  is  a  representation  of  the  velocity  structure  function  created  in  the  same 
manner.  By  making  the  assumption  that  the  flow  is  isotropic  and  homogeneous  there 
will  be  no  difference  in  this  function  no  matter  the  reference  orientation  is  made.  In 
Figure  6  it  is  important  to  note  behavior  at  the  extremes  of  the  position  scale.  As  the 
scale  size  becomes  small  there  is  very  little  variation  in  the  velocity  held  so  it  goes  to 
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zero  as  the  separation  gets  smaller.  As  the  separation  get  larger  the  function  tends  to 
become  constant  as  the  differences  in  the  spatial  parameter  start  to  average  out  to  a 
constant  value.  Using  an  r2^  fit  relationship  there  is  a  region  where  the  fit  coincides 
well  with  the  structure  function.  This  region  represents  the  inertial  subrange.  The 
structure  function  can  be  found  to  fit  within  this  region  by  the  following  relationship. 

A, Or)  =  C2vri  (10) 

In  Equation  10  two  important  things  are  presented.  The  first  is  the  r2/3  relationship. 
The  second  is  the  fit  constant  C2,  which  is  how  the  curve  fit  is  made.  This  Structure 
Constant  can  give  some  indication  of  how  relatively  weak  or  strong  the  fluctuations 
are. 


r(m) 

Figure  6.  Velocity  Structure  Function  and  Fit  Example.  [14] 


This  can  be  used  to  create  a  structure  function,  Dn(r ),  for  the  index  of  refraction 
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as  seen  below  in  Equation  11. 


Dn  (r ) 


J  C2nlo3\2,  0  <r</0 
[Cy*,  l0^r^:L0 


(11) 


The  index  of  refraction  structure  constant  (C‘l)  is  inferred  from  the  variability  of  the 
index  of  refraction  shown  in  Equation  11.  C2n  is  typically  in  the  range  of  10-12  for 
cases  of  strong  turbulence  to  10-17  in  weak  cases.  For  atmospheric  representations, 
Cl  is  typically  derived  from  C\  which  is  constructed  from  thermal  gradient  data 
collected  in  the  atmosphere  by  the  relationship  seen  in  Equation  12.  [2] 


C2n  =  C“( 79  x  10”6( 


_6  ,P(m&ars) 


T(K) 


)f 


(12) 


2.3. 2. 2  Fried  Parameter. 

In  trying  to  optimize  resolving  power  of  an  optical  system,  D.  L.  Fried  developed 
a  normalized  signal-to-noise  ratio  (^(^r))  as  a  function  of  the  aperture  diameter  (D) 
and  a  length  scale  (r0)  which  is  often  defined  as  the  Fried  Coherence  Length  or  Fried 
Parameter.  [5]  Figure  7  is  a  plot  of  this  function.  The  Fried  Parameter  is  found  at 
the  intersection  of  the  two  assymptotic  limits  at  ^  =  1.  The  Fried  Parameter  can  be 
found  from  the  structure  function  of  phase  and  amplitude  for  normalizing  the  Fried 
signal-to-noise  ratio  by  the  relationship  provide  below  in  Equation  13. 


D(j>{r )  +  Dln{A)(r)  =  6.88(— )s  (13) 

^  O 

The  Fried  Parameter  is  a  way  to  describe  the  point  of  diminishing  returns  for  improv¬ 
ing  signal  quality  by  increasing  the  aperture  diameter  to  create  a  diffraction  limited 
beam. 
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DIAMETER,  D/ro 

Figure  7.  Fried’s  Normalized  Signal-to-Noise  ratio.  [5] 

The  Fried  Parameter  can  be  found  with  a  given  using  the  expression  below  in 
Equation  14  for  a  plane  wave, 

2.905  2-IT  2  fR  2  3/5 

r°=^~6  88  ^T^J  C^s)ds\  (14) 

where  R  is  the  path  length  or  range.  If  C^(s)  does  not  vary  much  over  that  range 
then  it  is  possible  to  simplify  that  expression  to  Equation  15  for  a  plane  wave.  [11] 

As 

rQ  =  0.185  3  5  3  (15) 

R*C*  f 

2.4  Deviations  from  Kolmogorov 

By  now  understanding  the  origins  of  the  stochastic  representation  of  turbulence,  it 
is  important  to  note  that  there  are  several  problems  when  trying  to  capture  the  effects 
of  atmospheric  turbulence  on  beam  propagation  using  Kolmogorov  representation. 
The  first  is  the  fact  that  it  is  only  viable  in  the  inertial  subrange  of  eddy  structure. 
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The  second  is  that  it  assumes  a  homogeneous,  isotropic  distribution  of  turbulence. 
However  by  its  nature  turbulence  is  anything  but  homogeneous  and  isotropic.  Lo¬ 
calized  turbulent  disruptions  and  non-uniform  distributions  lead  to  anisotropy.  Tur¬ 
bulence  structures  are  not  only  circular,  but  have  helicity  to  them  which  can  form 
localized  patches  of  anisotropy. 

Golbraikh  et  al  found  that  using  Lidar  measurements  to  compare  vertical  turbu¬ 
lence  profiles,  that  Kolmogorov  theory  cannot  accurately  describe  turbulence  effects 
on  optical  beam  propagation.  [6]  Toselli  et  al  summarized  recent  experimental  results 
that  did  not  agree  with  analytical  techniques  utilizing  Kolmogorov  theory  for  atmo¬ 
spheric  turbulence.  Mechanisms  for  these  deviation  stem  from  density  anisotropy, 
helicity  of  turbulence,  and  atmospheric  boundary  layer  stability.  [18]  Shugaev  et  al 
found  that  the  Kolmogorov  theory  does  not  accurately  provide  the  correct  structure 
definitions  for  flow  fields  that  include  strong  disturbances  and  non-uniformity  of  the 
medium.  The  results  also  indicate  that  in  order  to  capture  turbulent  flow  accurately, 
it  is  necessary  to  have  a  large  number  of  data  points  in  order  to  get  the  appropriate 
structure  resolution.  [13]  Sirazetdinov  et  al  compared  the  results  from  an  experiment 
collecting  laser  beam  propagation  data  through  a  turbulent  engine  flow  field  with  that 
of  an  analytical  model  using  Fresnel  transformation  with  randomly  inhomogeneous 
phase  screens.  [15]  They  found  that  there  is  some  improved  accuracy  between  exper¬ 
iment  and  computations,  which  suggest  that  accounting  for  a  chaotic  turbulent  flow 
field  may  be  needed  to  accurately  capture  turbulence  effects  on  beam  propagation. 

2.5  Fluid  Studies  and  Aero-Optics 

There  have  been  several  studies  to  determine  the  behavior  of  turbulence  for  ap¬ 
plications  in  aero-optics.  Siegenthaler  et  al  provided  a  compilation  of  the  prevailing 
statistical  methods  and  parameters  with  the  assumptions  needed  to  apply  them  and 
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compared  them  to  the  effects  of  an  aero-optical  flow  held.  [14]  They  conclude  that 
making  stochastic  assumptions  are  not  viable  for  describing  aero-optical  hows  because 
Stochastic  methods  are  primarily  driven  by  thermal  gradients  while  aero-optical  tur¬ 
bulence  is  driven  primarily  by  how  gradients  in  the  shear  layer.  They  do  say  that  these 
statistically  derived  parameters  could  be  viable  on  a  case  by  case  basis  for  detailed 
scenarios.  Zubair  et  al  provides  two  important  perspectives  to  this  research  topic. 
[22]  The  hrst  is  that  the  large  scales  of  turbulent  structure  provide  more  dominant  ef¬ 
fect  on  aero-optical  performance.  The  second  is  that  a  Large  Eddy  Simulation(LES) 
turbulence  model  provides  enough  resolution  to  evaluate  aero-optical  performance. 
High  Reynolds  number  how  helds  are  dominated  by  smaller  scale  turbulent  eddies. 
While  large  scale  eddies  have  more  deleterious  effects  on  beam  propagation,  small 
eddies  contribute  to  beam  degradation  as  well,  but  are  slightly  more  difficult  to  mea¬ 
sure  and  predict.  Catrakis  et  al  took  the  approach  of  combining  CFD  research  with 
an  experimental  method  to  observe  fluid  turbulence.  [3]  The  work  compared  exper¬ 
imental  results  with  computational  simulations  to  observe  aero-optical  effects.  The 
major  folding  was  that  the  turbulence  behavior  is  non-monotonic,  which  means  that 
there  are  non-linear  effects  in  turbulence  that  may  not  be  able  to  be  resolved  purely 
statistically.  Aguirre  et  al  found  that  at  higher  Reynolds  numbers,  the  number  of 
small-scale  features  increases  which  creates  the  inability  to  compare  the  propagation 
effects  from  large  scale  with  regard  to  small  scale  structures.  Small  scales  are  consis¬ 
tent  with  similarity  theory  for  statistical  means  of  representation.  [1]  However  there 
still  exists  a  large  anisotropy  of  the  density  which  must  be  compensated  for  when 
computing  wave  front  propagation. 

Visbal  researched  effects  on  optical  propagation  through  the  shear  layer  of  turbu¬ 
lence  for  various  airspeeds  using  two  dimensional  high-order  numerical  fluid  methods. 
[19]  The  conclusions  found  that  transmission  aberrations  compared  using  Optical 
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Path  Difference  can  manifest  cine  to  primarily  compressibility  effects,  despite  very 
different  density  field  distributions.  This  suggests  that  there  is  no  one-to-one  cor¬ 
respondence  between  integrated  optical  aberration  and  the  flow  structure,  meaning 
that  there  are  several  mechanisms  contributing  to  adverse  propagation  effects  through 
turbulence. 

These  studies  provide  some  interesting  insight  into  the  shortcomings  of  Kol¬ 
mogorov  theory  and  the  difference  in  approaches  between  atmospheric  and  aero-optics 
turbulence  communities.  Kolmogorov  is  useful  for  atmospheric  turbulence  in  that 
provides  a  robust  method  for  accessing  the  affects  of  the  atmosphere  on  beam  prop¬ 
agation  despite  slight  deviation  from  the  prevailing  theory.  Atmospheric  turbulence 
community  also  uses  thermal  gradients  in  their  determinations  as  it  is  the  primary 
source  for  atmospheric  turbulence  and  is  data  that  is  easily  obtained  from  rnetero- 
logical  instrumentation.  Aero-optical  turbulence  consists  of  turbulent  scales  that  are 
smaller  in  scale  driven  by  flow  gradients  in  the  shear  layer.  These  eddies  are  suscepti¬ 
ble  to  compressibility  factors  creating  density  gradients  changing  the  refractive  index 
as  opposed  to  thermal  gradients.  The  largest  scales  of  turbulence  in  aero-optics  flows 
are  the  most  dominant  influence  on  beam  propagation,  and  the  ability  to  manage  or 
decrease  their  size  will  have  less  critical  effects  on  beam  propagation.  This  is  relevant 
to  this  research  by  investigating  if  these  differences  in  approach  are  quant  ifiably  dif¬ 
ferent,  and  that  even  if  they  are,  will  the  difference  significantly  alter  what  method 
is  chosen  to  analyze  an  aero-optic  flow’s  affect  on  beam  propagation. 

2.6  CFD  Tools  and  Use  in  Aero-Optical  Investigations 

Modeling  numerical  fluid  flow  can  be  accomplished  using  the  Navier-Stokes  equa¬ 
tions.  Numerical  modeling  of  fluid  dynamics  using  the  Navier-Stokes  equations  has 
become  very  robust,  and  there  are  many  commercially  available  flow  solvers  that 
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allow  these  flow  physics  to  be  calculated.  As  understanding  of  turbulence  behavior 
has  evolved,  mathematical  models  that  predict  the  unsteady  nature  of  turbulence 
have  been  incorporated  with  these  Navier-Stokes  solvers  to  provide  a  characteristic 
representation  of  the  random  nature  of  turbulence. 

One  of  these  is  the  Government  developed  code,  the  Air  Vehicles  Unstructured 
Solver  (AVUS).  AVUS  is  an  unstructured,  cell-centered,  finite- volume,  Godunov-type 
solver  that  uses  least-squares  gradient  reconstruction  and  limiting  for  second-order 
spatial  accuracy,  and  second-order,  point-implicit  time  integration.  [7;  16]  The  com- 
mericial  equivalent  of  AVUS  is  Cobalt.  Rennie  et  al  used  Cobalt  to  compare  a  2-D 
compressible  shear  layer  to  a  Weakly  Compressible  Model.  [12] 
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III.  Methodology 


This  chapter  presents  the  reader  with  the  steps  taken  to  conduct  this  investigation. 
The  primary  scope  of  this  research  is  to  develop  a  methodology  by  which  to  make  these 
assessments  with  a  variety  of  scenarios.  As  an  initial  test  case  the  flow  scenario  will  be 
a  2-D  laminar  Navier-Stokes  solver  weakly  compressible  aero-optical  like  flow  field.  In 
order  to  create  the  basis  for  this  comparison  of  turbulence  effects,  a  fluid  flow  similar 
to  an  aero-optic  flow  will  be  generated  using  a  laminar  Navier-Stokes  solver  used  in 
CFD  software.  A  2-D  splitter  plate  type  geometry  is  used  to  generate  a  laminar  shear 
layer  as  representative  turbulent  field  for  analysis.  Stochastic  methodologies  will  then 
be  used  to  assess  this  output  and  determine  if  this  characteristic  aero-optical  turbulent 
flow  field  has  Kolmogorov  like  qualities. 

Figure  8  provides  a  visual  representation  of  the  workflow  for  this  research.  There 
are  three  primary  phases  to  accomplishing  this.  The  first  is  pre-processing  or  the 
preparing  of  inputs  for  processing.  A  geometry  configuration  is  selected  and  a  carte¬ 
sian  grid  is  constructed.  Combining  this  grid  with  the  boundary  conditions  and  the 
run  script  the  CFD  process  may  begin  which  is  one  of  the  components  for  the  sec¬ 
ond  phase,  Processing.  Processing  is  the  computation  of  the  flow  and  calculation  of 
the  optical  parameters  from  the  flow  data.  The  last  phase  is  Post-Processing.  Post- 
Processing  is  the  required  steps  to  examine  the  flowficld  and  assess  that  the  data  is 
either  complete  or  sound,  and  then  to  analyze  the  statistical  parameterization. 

3.1  Geometry 

The  geometry  to  be  used  for  this  investigation  is  a  2-D  splitter  plate.  Figure  9 
is  an  exaggerated  illustration  of  this  geometry.  A  splitter  plate  is  a  thin  sheet.  The 
flow  above  and  below  are  set  to  varied  speeds.  When  the  flow  reaches  the  end  of  the 
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Figure  8.  Work  Flow  Chart  for  the  Investigation 

plate  the  shear  induced  turbulence  will  evolve.  For  this  case,  the  freestream  Mach 
number  (M^)  above  the  plate  is  set  to  0.6  (Mij00  =  0.6).  Below  the  plate  it  is  set 
to  0.1  (M2j oo  =  0.1).  This  is  consistent  with  making  the  laminar  shear  layer  weakly 
compressible  and  is  similiar  to  the  values  set  in  Visbal  and  Rennie  et  al.  [19;  12] 

3.2  Grid 

A  cartesian  grid  was  created  using  the  software  package  Gridgen.  A  structured  grid 
was  used  due  to  the  simplicity  of  this  geometry,  improved  control  on  mesh  spacing, 
and  the  fact  that  my  statistical  parameterization  algoritm  uses  a  structured  grid. 

Before  building  the  grid  it  was  necessary  to  make  some  calculations  regarding  the 
minimum  grid  spacing  at  the  end  of  the  plate  so  as  to  be  able  to  capture  the  flow 
physics  within  the  boundary  layer.  The  boundary  layer  thickness  (5)  at  the  end  of 
the  plate  can  be  found  by  using  Equation  16  for  a  laminar  boundary  layer  as  seen 
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Figure  9.  Work  Flow  Chart  for  the  Investigation 

below. 

0.37/p 
^  1 
Re  3 

For  the  purposes  of  defining  <5,  the  Reynold’s  Number  was  set  to  4000  (Re  =  4000). 
This  Reynold’s  Number  was  used  because  the  flow  solver  is  using  the  Laminar  Navier- 
Stokes  equantion  and  a  Reynold’s  Number  of  4000  is  consistent  with  previous  research 
for  the  purpose  of  a  laminar  flow  analysis.  In  defining  the  geometry  the  splitter  plate 
length  (lp)  is  set  to  lm  to  better  allow  the  flow  to  develop.  From  these  parameters 
the  boundary  layer  thickness  is  found  to  be,  6  =  0.023m.  This  was  determined  by 
using  the  top  layer  free  stream  velocity  and  assuming  that  the  lower  layer  free  stream 
velocity  would  generate  a  similiar  boundary  layer  thickness  as  was  used  in  previous 
research. 

In  order  to  include  10  grid  points  with  in  the  boundary  layer  using  a  geometric 
growth  rate  of  1.2  for  the  node  spacing  a  minimum  grid  spacing  (d')  was  found  to  be 
0.00075m.  The  width  of  the  splitter  plate  was  set  to  26'  in  order  to  the  plate  thin. 
From  these  calculations  the  grid  was  created.  Figure  10  is  a  view  of  the  whole  grid. 
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The  entire  grid  has  a  dimension  of  6 m  x  7m.  The  grid  dimensions  are  exagerated  so  as 
to  try  to  eliminate  interference  within  the  area  of  interest  from  reflections  and  other 
boundary  influences.  The  detail  within  the  area  is  sufficient  enough  so  as  to  capture 
the  flow  physics.  Figure  11  shows  the  grid  zoomed  in  to  the  edge  of  the  splitter  plate. 
This  view  gives  an  idea  of  the  growth  and  mesh  size  based  on  the  above  definitions. 


Figure  11.  Zoomed  in  View  of  the  Grid 
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3.3  Boundary  Conditions 


The  boundary  conditions  tell  the  flow  solver  what  conditions  are  to  be  held  con¬ 
stant  or  parameterize  the  behavior  of  the  flow  along  the  geometry  boundaries.  The 
inlet  regions  were  set  to  source  boundaries  as  they  will  be  where  the  flow  will  be  in¬ 
troduced  from.  The  parameters  of  these  boundaries  are  built  on  assuming  Standard 
Temperature  and  Pressure  (STP)  with  the  Mach  numbers  for  the  upper  and  lower 
surfaces  set  to  values  specified  in  the  geometry  session.  The  top  and  bottom  bound¬ 
aries  are  set  to  slip  walls  so  as  to  bound  the  flow  but  not  create  another  boundary 
layer  that  would  be  introduced  if  it  was  a  solid  wall.  The  outflow  region  is  set  to 
a  farfield  boundary  as  it  is  sufficiently  downstream  from  the  area  of  interest  so  as 
not  to  affect  the  flow  there.  The  splitter  plate  was  set  to  an  isothermal  no-slip  wall. 
Isothermal  was  chosen  so  as  to  not  influence  the  boundary  layer  development  through 
thermal  fluctuations.  The  boundary  condition  file  used  for  this  flow  scenario  can  be 
found  in  Appendix  A. 

3.4  Job  File 

The  job  file  is  the  scripting  that  initializes  the  flow  solver.  Within  it  are  the 
commands  which  allow  for  running  it  parallel  and  an  input  file  containing  the  solver 
parameters.  The  major  things  to  note  here  are  the  run  conditions  and  the  reference 
parameters.  The  solver  run  conditions  were  set  to  Laminar  Navier-Stokes  with  2000 
iterations.  Time  accuracy  of  the  solver  was  turned  off  and  a  large  CFL(CFL=106) 
was  used  to  have  the  time  scale  adjust  to  allow  the  flow  to  develop  quickly  to  a 
representative  field.  2000  iterations  were  chosen  so  as  to  allow  this  flow  to  converge 
to  representative  turbulent  field.  The  resultant  flow  field  represents  a  snapshot  of 
what  a  numerically  resolved  field  should  resemble,  but  has  not  used  specified  local 
time-stepping  for  continual  flow  evolution.  The  spatial  resolution  was  set  to  second 
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Table  1.  Flow  Reference  Parameters 


Parameter 

Value 

Refrence  Pressure  (P0) 

P0  =  101325 Pa 

Refrence  Temperature  (T0) 

T0  =  288.15 K 

Refrence  Top  Inlet  Mach  Number  (M1>00) 

Tfjoc  =  0.6 

Refrence  Bottom  Inlet  Mach  Number  (M2i00) 

M2)00  =  0.1 

order  while  the  temporal  resolution  was  set  to  first  order  because  the  purpose  is  to 
get  a  steady-state  flow  field.  The  job  hie  used  for  this  how  scenario  can  be  found  in 
Appendix  B. 

3.5  Fluid  Solver 

AVUS  was  used  for  the  CFD  work  for  this  scenario.  AVUS  is  an  Eulcr/Navier- 
Stokes  solver,  was  used  to  obtain  solutions  on  these  grids.  AVUS  is  an  unstructured, 
cell-centered,  finite-volume,  Godunov-type  solver  that  uses  least-squares  gradient  re¬ 
construction  and  limiting  for  second-order  spatial  accuracy,  and  second-order,  point- 
implicit  time  integration.  [7;  16]  It  handles  two  and  three  dimensions,  arbitrary  cell 
types,  and  has  been  efficiently  parallelized  using  Message  Passing  Interface  (MPI). 
AVUS  was  used  because  of  previous  successful  work  in  this  type  of  research,  and  the 
author’s  familiarity  with  the  program.  For  AVUS  to  work  it  must  hrst  be  compiled 
for  the  parallel  processing  enviornment  of  the  computational  resources  being  used. 
AVUS  was  compiled  on  the  Linux  clusters  in-house  at  AFIT.  The  job  hie  contains  all 
the  scripting  required  to  run  AVUS  on  the  AFIT  clusters. 

3.6  Fluid  Computation  Results  Interpretation 

At  the  completion  of  the  how  solver,  an  intermediate  post-processing  step  is  com¬ 
pleted  to  determine  the  completion  and  convergence  to  a  completed  representative 
how  held  that  is  steady  state.  The  hrst  step  is  to  examine  the  completed  how  held 
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using  a  flow  visualization  package.  Fieldview  is  an  example  of  a  software  which  can 
do  this.  The  grid,  boundary  conditions,  and  the  resultant  data  are  loaded  into  the 
software.  This  is  useful  in  that  it  provides  a  visual  diagnostic  tool  to  examine  the  flow 
held  for  any  errors  or  problematic  behavior.  If  there  are  problems  it  may  be  neces¬ 
sary  to  adjust  the  solver  inputs  or  grid/boundary  hies  to  correct  for  these  problems. 


Figure  12.  Computed  Flow  Field  Vorticity  Magnitude 


Figure  13.  Zoomed  in  View  of  the  Vorticity  Magnitude  at  the  end  of  the  Splitter  Plate 
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Figure  12  provides  a  view  of  the  converged  flow  field  vorticity  magnitude  used  for  this 
flow  scenario.  In  Figure  12  there  is  some  noise  present  in  the  downstream  boundary 
and  what  looks  like  a  bubble  in  the  lower  half,  midway  through  the  flow  held.  This 
could  be  sources  introduced  by  the  boundary  or  grid  resolution.  However,  these  ar¬ 
tifacts  were  ignored  because  it  was  outside  of  the  region  of  interest  for  the  analysis. 
Figure  13  gives  a  more  in  depth  look  at  the  vorticity  at  the  end  of  the  splitter  plate 
(the  area  of  interest).  It  is  evident  to  note  that  there  is  very  strong  circulation  at  the 
end  of  the  plate  which  is  expected.  It  is  also  evident  that  the  resolution  seems  to  be 
good  enough  to  be  capturing  the  eddy  structures  coming  off  the  plate. 

Convergence  is  determined  by  analyzing  the  trends  in  iteration  residuals  which 
AVUS  generates  as  a  secondary  output.  Figure  14  shows  the  residual  convergence 
trend  for  the  test  case.  The  residual  started  to  level  off  at  about  1500  iterations.  It 
appears  to  still  be  decreasing,  however  the  trend  suggest  that  it  should  be  converged 
to  point  where  further  solver  iterations  will  provide  only  marginal  improvements  to 
the  flow  field  solution  data.  The  Y+  is  a  non-dimensional  representation  of  the  wall 
spacing.  For  turbulent  cases  Y+  needs  to  be  as  close  to  one  as  possible  depending 
on  the  boundary  conditions.  Figure  15  is  the  convergence  trend  for  Y+  for  this  case 
as  generated  by  AVUS.  It  appears  to  be  settling  around  a  Y+  in  the  low  forties. 
This  means  that  this  grid  will  need  to  be  resized  to  be  used  in  generating  turbulent 
Navier-Stokes  solutions.  Here  Y+  is  used  as  a  diagnostic  tool  to  check  for  convergence. 

3.7  Exporting  Flow  Data 

Following  the  determination  of  a  converged  steady-state  laminar  shear  layer  flow 
field,  it  is  important  to  get  the  actual  data  values  calculated  from  the  CFD  solver. 
The  Blacksmith  utility  code  developed  for  use  with  AVUS  grids  and  results  has  the 
capability  of  generating  this  geometric  data  and  their  corresponding  flow  parameter 
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values  in  a  formatted  tabular  file.  It  is  this  file  that  is  imported  into  the  developed 
optical  parameterization  algorithm. 


Iterations 


Figure  14.  CFD  Solution  Residual  Convergence  Trend 


Iterations 


Figure  15.  CFD  Solution  Y+  Convergence  Trend 
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3.8  Optical  Parameterization  Algorithm 


A  Matlab  code  was  developed  for  calculating  the  parameters  needed  for  the  anal¬ 
ysis  for  this  study.  The  algorithm  reads  in  the  flow  data  and  generates  a  focused 
mesh  of  this  data  for  the  area  of  interest.  It  then  calculates  numerous  optical  pa¬ 
rameters  discussed  in  Chapter  2  for  the  development  of  structure  functions  and  the 
corresponding  structure  constants.  This  algorithm  can  be  found  in  Appendix  C. 

It  is  important  to  highlight  here  some  of  the  defining  values,  area  of  interest,  and 
path  defintions.  Table  2  provides  a  list  of  some  of  the  optical  parameter  values  used  in 
the  approach.  The  region  this  analysis  was  conducted  on  was  a  lm  x  1  m  box  centered 
on  the  end  of  the  splitter  plate  (A"  =  0  is  the  end  of  the  splitter  plate).  A  uniform 
mesh  was  constructed  in  this  box  using  a  spacing  of  0.001  giving  1001  x  1001  data 
points  for  the  calculations.  The  reference  point  for  the  calculation  of  the  structures 
was  set  to  the  spatial  values  set  at  X  =  0.5 m  and  Y  =  0 m  which  is  the  center  of  the 
mesh.  Because  the  structure  functions  are  calculated  from  spherical  coordinates  it  was 
necessary  to  calculate  the  radial  distances  from  the  mesh  center.  A  tolerance  was  used 
in  the  radial  distance  determination  for  calculating  the  structure  functions.  A  plane 
wave  was  assumed  for  the  beam  as  this  is  a  2-D  flow  held.  For  the  path  dependent 
integrals  and  aperture  averaging  the  Y-direction  represents  the  propagation  direction 
and  the  X-direction  represents  the  aperture  face.  A  visualization  of  this  approach  can 
be  seen  in  Figure  16. 

3.9  Output 

The  data  is  calculated  in  Matlab.  Matlab  has  great  plotting  capabilities  and  was 
used  to  generate  data  plots  for  interpretation.  A  more  in-depth  presentation  of  this 
can  be  found  in  Chapter  4. 


35 


Figure  16.  Parameterization  Algorith  Assumptions  Visual 


Table  2.  Reference  Data  for  Parameterization  Algorithm 


Beam  Waist  (w0) 

wQ  =  0.5m 

Amplitude  Peak  (A0) 

Aa  =  l 

Wavelength  (A) 

A  =  lm 

Reference  Density  (p0) 

Po  =  1.2215 

Gladstone-Dale  Constant  ( KGD ) 

KGd  =  2.25  x  10-4^ 

Radius  Tolerance 

tol  =  0.0005m 

Grid  Spacing  (Ax,  Ay) 

Ax  =  0.001m,  Ay  =  0.001m 

Radial  Stepsize  (Ar) 

Ar  =  0.001m 
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IV.  Results  and  Analysis 


Chapter  4  is  a  compilation  of  the  results  obtained  for  defining  the  stochastic 

optical  paramters  for  this  2-D  laminar  shear  layer  for  the  splitter  plate  geometry.  The 

results  begin  with  an  analysis  of  the  flow  field,  with  regard  to  the  actual  Reynold’s 

Number  within  the  mixing  layer  and  require  scale  resolution.  This  leads  to  the  the 

development  of  the  local  refractive  indicies  and  OPD’s  with  a  comparison  between 

the  two  methods  used  for  calculating  the  refractive  index  in  these  developments.  This 

is  followed  by  the  presentation  of  the  results  of  the  structure  function  development 

2 

and  the  constants  calculated  assuming  the  r  3  power  law  including  a  discussion  on 
deviations  from  expectations  and  reasonableness  of  the  numbers  calculated. 


Figure  17.  Area  of  Interest  Density  Mesh 
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Figure  18.  Streamwise  Reynold’s  Number  within  the  Mixing  Layer  for  the  Area  of 
Interest 


-5 


Figure  19.  Streamwise  Smallest  Scale  Size  Resolution  within  the  Mixing  Layer  for  the 
Area  of  Interest 

4.1  Flow  Field  in  Area  of  Interest 


As  described  in  Chapter  3  the  area  of  interest  is  confined  to  1  m  x  1  m  box  centered 
at  the  end  of  the  splitter  plate.  Density  is  the  primary  flow  parameter  in  determining 


the  local  refractive  indicies.  Figure  17  is  a  view  of  the  local  densities  within  this  area 
of  interest  as  meshed  by  Matlab.  This  figure  gives  a  good  view  of  the  eddies  flowing 
from  the  end  of  the  plate  with  regard  to  variations  from  higher  density  (red)  to  lower 
density  (blue). 

The  flow  field  seems  to  capture  some  of  the  large  scale  features,  but  doesn’t  really 
provide  any  insight  into  if  the  flow  in  the  mixing  layer  is  turbulent  or  not.  The  best 
way  to  investigate  this  is  to  calculate  the  Reynold’s  Number  along  the  streamwise 
direction  within  this  area  of  interest.  The  Reynold’s  Number  is  calculated  using 
Equation  1.  The  flow  speed  for  this  analysis  is  calculated  by  the  following  equation 
where  the  Mixing  Layer  Flow  Speed(AfZ)  is  the  average  of  the  upper  freestream 
velocity  (£4)  and  the  lower  free  stream  velocity  (£//). 

A U  =  V h  ~  ^  (17) 


The  length  scale(L(x))  used  is  the  half  width  of  the  mixing  layer  which  is  found  by 
taking  the  average  in  difference  between  the  Mixing  Layer  Flow  Speed  (A U)  and  the 
Lower  Freestream  Velocity  (Ui). 


L(x) 


A  U-Ut 
2 


(18) 


The  streamwise  Reynold’s  Number  can  be  found  in  Figure  18.  The  Reynold’s  Number 
grows  linear  within  this  area  of  interest  and  quickly  approaches  105  which  shows  that 
this  mixing  layer  is  turbulent  and  should  have  a  Kolmogorov  like  distribution  of  scales. 

However,  the  grid  was  constructed  assuming  a  laminar  shear  layer  so  it  is  necessary 
to  compute  what  the  smallest  resolution  is  required  to  capture  the  smallest  eddy 
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scales.  The  inner  scale  (l(x))  is  computed  by  using  the  following  relation. 


l(x) 


(0.1  Rel) 


(19) 


This  result  is  plotted  in  Figure  19.  This  shows  that  the  smallest  scale  is  on  the  order 
of  10_5m.  However  the  smallest  scale  size  for  this  grid  is  10~3m.  This  means  that 
while  the  mixing  layer  is  turbulent  with  regard  to  Reynold’s  Number  the  scale  sized 
used  may  not  effectively  capture  the  entire  distribution  of  scales,  which  may  limit  the 
feasibility  in  the  interpretation  of  the  results.  However,  there  are  sufficient  large  scale 
eddies  present  which  have  a  more  dominant  effect  on  optical  beam  propagation. 


4.2  Refractive  Index  and  OPD 

Using  the  density  mesh  presented  above  the  local  refractive  indicies  were  calculated 
and  meshed  using  the  two  methods  described  in  Chapter  2:  Edlen  (Equation  2)  which 
includes  wavelength  and  Gladstone-Dale  (Equation  3).  Figures  20  and  21  show  the 
mesh  calculated  using  the  above  equations.  A  quick  comparison  of  the  two  figures 
reveals  that  the  two  fields  are  indistinquishable.  This  suggests  that  at  1  micron 
wavelength  the  calculations  are  very  close.  The  other  observation  that  can  be  made 
is  that  the  index  variation  maps  very  closely  to  the  eddies  in  the  flow  as  expected 
since  they  are  density  dependent  and  are  evident  in  Figure  17. 

The  local  OPD  is  calculated  from  the  local  OPL’s  generated  by  the  local  refractive 
indicies  as  described  in  Equation  6.  Figures  22  and  23  show  the  OPD  mesh  for  the 
respective  fields.  Much  like  the  refractive  index  meshes  these  are  indistinquishable 
which  is  consistent  with  the  fact  that  the  calculations  are  based  on  the  refractive 
index.  Figure  24  provides  a  sample  of  the  OPD  at  the  top  of  the  box  along  the 
streamwise  direction  for  the  two  different  OPD’s.  They  seem  to  match  very  well. 
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Figure  20.  Local  Refractive  Index  Using  Edlen  Formula 


Figure  21.  Local  Refractive  Index  using  Gladstone-Dale 
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Figure  22.  Local  OPD  from  Edlen  Refractive  Indicies 


Figure  23.  Local  OPD  from  Gladstone-Dale  Refractive  Indicies 
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Figure  24.  OPD  Comparison  between  Refractive  Index  Methods 


X  10A-8 


Figure  25.  OPD  Absolute  Difference  of  Methods 


Figure  25  provides  an  absolute  difference  in  the  two  values.  The  difference  between 
them  on  the  order  of  8  decimal  places  which  suggests  there  is  negligible  difference  in 
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the  two  refractive  index  methods  for  this  wavelength.  This  provides  a  reasonable 
argument  for  choosing  a  single  method  for  these  calculations  at  this  wavelength. 
However,  both  methods  were  used  throughout  this  algorithm  so  as  to  ensure  similarity 
of  methods. 

4.3  Structure  Functions  and  Constants 

The  first  structure  function  calcuted  was  Dn(r )  using  the  method  described  in 
Equation  9  directly  from  the  flow  field  data.  The  results  from  this  calculation  can 
be  found  in  Figure  26.  The  data  trend  should  resemble  Figure  6  in  that  it  should 
go  rapidly  to  zero  as  the  separation  distance  (r)  gets  smaller.  This  is  not  what  is 
seen  here.  This  could  be  a  result  of  a  lack  of  datapoints  for  sampling  at  the  smallest 
scales.  Figure  27  gives  the  number  of  datapoints  sampled  for  all  the  separation 
distances  for  the  averages.  This  gives  some  indication  of  the  resolution  limiting  at 
the  smallest  scales  which  prevents  the  trend  from  happening  as  expected.  At  the 
largest  scales  there  appears  to  be  some  oscillation  in  the  structure  function.  These 
could  be  due  to  the  selection  of  the  reference  point  in  that  there  may  be  stronger 
variations  in  the  refractive  index  respective  to  that  reference  point.  Because  this  was 
a  single  timestep  and  a  single  reference  point  there  could  be  a  lack  of  sampling  data 
in  these  regions  inducing  these  oscillations.  This  might  be  able  to  be  corrected  using 
various  references  and  averaging  or  using  time  averaged  data  for  a  time  evolving  flow. 
This  also  demonstrates  that  this  flow  is  anything  but  isotropic  and  homogeneous. 
However  these  limitation  could  be  due  to  the  coarseness  in  the  grid  resolution  which 
can’t  capture  the  smaller  scale  eddies. 

From  the  refractive  index  structure  function,  Dn(r),  the  C ®  was  calculated  using 
2 

the  r  s  power  law.  This  can  be  seen  in  Figure  28.  These  points  are  the  circles(Edlen) 
and  x(Gladstone-Dale).  These  C„  values  represent  what  they  would  be  if  Kolmogorov 
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Figure  26.  Refractive  Index  Structure  Function  as  Calculated  from  Varying  Refractive 
Index’s  in  the  Flow  Field 


Figure  27.  Number  of  Data  Points  Sampled  for  Corresponding  Separation  Distance 
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2 

was  assumed  using  the  r  3  power  law.  The  first  thing  to  note  is  that  this  not  a  constant 
2 

using  the  r  3  power  law  which  suggests  that  the  aero-optical  flow  is  non- Kolmogorov. 
The  other  thing  is  the  range  of  C2’ s  is  between  1CT12  and  1CT10  which  is  higher  than 
expected  for  definitions  of  atmospheric  turbulence  (Typical  Range:  1CT17  <  < 

ICE12).  This  seems  to  suggest  that  this  particular  flow  field  is  non- Kolmogorov. 
Again  though  this  was  a  single  data  field  for  a  specific  case.  Time  averaged  data 
or  varied  reference  points,  a  more  finely  resolved  grid,  or  a  turbulent  Navier-Stokes 
produced  case  may  yield  more  Kolmogorov  like  turbulent  profiles. 
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Figure  28.  Calculated  Structure  Constant  (C£)  from  Dra(r)  using  ry  Cf.  from  Flow 
Field,  and  rD  from  the  Flow  Field 

Figure  28  also  compares  various  other  methods  for  calculating  the  C2  from  other 
parameters.  The  red  +  signs  represent  the  C*2  values  found  by  assuming  Kolmogorov 
distribution  in  the  flow  field  for  the  local  temperatures.  Cj.  was  found  by  construct¬ 
ing  the  temperature  structure  function  (Dt{t)).  This  is  then  used  to  calculate  C2 
using  Equation  12.  The  asterisks  and  diamonds  is  the  C*2  as  calculated  using  Equa¬ 
tion  15  from  which  the  ra  values  were  calculated  from  the  field  using  Equation  13. 
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What  is  interesting  to  note  is  the  regions  where  the  trend  is  coincidental.  Building 
Cl  from  C %  overlaps  in  the  range  of  0.01m  —  0.1m  which  is  within  the  theoretical 
atmospheric  inertial  subrange.  It  is  also  interesting  to  note  that  using  thermal  gra¬ 
dients  for  Kolmogorov  Theory  corresponds  with  the  flow  field  derived  parameters 
over  a  longer  range  in  that  atmospheric  analysis  use  thermal  gradients  to  determine 
the  parameters.  The  agreement  by  calculating  C\  from  r0  has  a  much  smaller  range 
(0.01m  —  0.02m),  yet  is  also  within  the  theoretical  atmospheric  inertial  subrange. 

Figure  29  presents  the  computation  of  the  Fried  Parameter  for  two  methods.  The 

first  two  lines(Edlen  and  Gladstone-Dale)  represent  constructing  the  Fried  Parameter 

based  on  the  Phase  and  Amplitude  structure  functions  as  found  using  Equation  13. 

The  second  two  lines  (Cl)  represent  the  calculation  from  C2n  values  found  using  the 
2 

rs  power  law  with  Dn(r).  Again  the  values  assume  a  variable  range  rather  than  a 
dehnate  value.  However  from  a  separation  distance  range  of  0.01m  —  0.02m  there 
seems  to  be  some  agreement  for  rQ  ~  0.05m  which  is  a  reasonable  value  for  the  Fried 
Parameter.  It  is  also  interesting  to  note  the  differences  between  the  two  methods. 
This  can  be  seen  in  Figure  30.  There  is  a  range  of  values  from  which  they  differ 
however  there  is  an  average  difference  of  about  5 cm.  That  means  that  for  sizing 
r0  there  is  almost  insignificant  difference  in  value  between  a  statistical  representation 
and  the  actual  flow  field  derived  values  for  the  purposes  of  system  design  for  Adaptive 
Optics. 

Another  method  of  comparison  of  the  results  is  the  corelation  of  the  phase  struc¬ 
ture  function  (D^r))  calculated  from  the  field  and  the  D^r)  from  the  C2n .  These 
results  can  be  seen  in  Figure  31.  The  Edlen  and  Gladstone-Dale  points  represent  the 
values  calculated  from  the  flow  held.  The  Cl  lines  are  derived  from  using  the  values 
calculated  using  Kolmogorov.  From  this  figure  the  discussion  relative  to  the  function 
trend  presented  above  for  Dn(r)  applies  here. 
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Figure  29.  Calculated  Fried  Parameter  (rc)  from  D,f,,Din^  and  C2 
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Figure  30.  Calculated  Fried  Parameter  ( rQ )  Difference  between  Methods 
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Figure  31.  Phase  Structure  Function  Comparison  (D^,(r))  as  found  from  Flow  Field  and 
Kolmogorov  Derived  C 2 


Tabic  3  is  a  synopsis  of  the  interesting  details  from  these  results.  This  0.01m  — 
0.02m  range  is  very  interesting  in  its  repetition  and  the  fact  that  it  is  within  the 
atmsopheric  inertial  subrange.  This  may  be  purely  coincidental  though  without  fur¬ 
ther  grid  refinement  to  ensure  capture  of  the  smallest  turbulence  scales.  While  there 
are  major  quantifiable  differences  between  the  methods  for  this  shear  mixing  layer, 
this  suggests  that  further  scenarios  would  need  to  be  conducted  to  draw  any  relevant 
conclusions  on  if  this  is  purely  coincidental  or  has  some  applicable  ramifications. 


Table  3.  Summary  of  Interesting  Parameterization  Data 


Comparison 

Agreement  Value/Range 

rQ 

r0  «  0.05m 

r2 

io-12  -  -10“1U 

C%(r) from  Flow  to  C%(r)  from  C f. 

0.01m  <  r  <  0.1m 

C%(r) from  Flow  to  C%(r)  from  ra 

0.01m  <  r  <  0.02m 

rQ from  Flow  to  rQ  from 

0.01m  <  r  <  0.02m 

49 


V.  Conclusions 


This  chapter  begins  by  providing  a  cursory  summary  of  the  results.  It  then  out¬ 
lines  opportunities  for  improvement  to  the  existing  methodology.  It  concludes  with 
suggestions  regarding  future  work  in  the  area  of  research. 

5.1  Results  Summary 

Siegenthaler  et  al  concluded  that  aero-optic  flows  are  non-Kolmogorov  like  because 
they  suggest  that  centrifigal  forces  associated  with  the  eddies  are  more  dominant  than 
the  thermal  variations  in  these  eddies.  Stochastic  methods  work  in  the  atmosphere 
because  temperature  gradients  are  the  primary  force  behind  inertial  subrange  turbu¬ 
lence  paired  with  high  Reynolds  numbers  to  make  Kolmogorov  theory  applicable.  In 
aero-optical  flows  there  are  temperature  fluctuations  in  the  eddies,  but  the  primary 
source  of  turbulent  kinetic  energy  is  in  the  inertia  of  the  whirling  flow,  where  the 
whirls  produce  density  and  refractive  index  gradients  via  centrifigual  effects  which 
seems  to  be  consistant  with  what  Siegenthaler  et  al  were  suggesting.  Figures  28 
shows  that  the  Kolmogorov  derived  and  C\  from  C\  have  a  larger  range  of  over¬ 
lap  which  suggests  that  temperature  gradients  are  the  dominant  parameter  in  scaling 
for  turbulence  strength  using  Kolmogorov.  When  calculating  C\  from  a  flow  field 
characteristic  rQ  there  is  much  less  agreement  which  suggests  that  thermal  gradients 
were  not  the  dominant  parameter.  It  doesn’t  make  a  case  for  what  parameter,  but 
Siegenthaler  et  al  hypothesize  centrifugal  force  seems  to  be  a  relevant  factor.  In  light 
of  these  discussions,  the  2-D  shear  mixing  layer  aero-optic  like  flow  was  not  Kol¬ 
mogorov  despite  the  fact  that  the  mixing  layer  appears  to  be  turbulent.  This  could 
be  a  result  of  not  having  enough  grid  resolution  to  capture  the  smallest  eddy  scales. 

However,  stochastic  methods  might  serve  as  a  useful  way  to  derive  system  design 
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parameters  for  Laser  Weapons  Systems.  The  Fried  Parameter  is  often  used  to  deter¬ 
mine  the  acuator  numbers  and  scaling  for  the  Adaptive  Optics  system.  By  comparing 
Kolmogorov  based  analysis  with  that  of  flow  field  measurement  there  were  differences, 
but  these  differences  were  small  (approximately  5 cm).  While  this  isn’t  insignificant,  it 
suggests  that  a  stochastic  method  provides  a  reasonable  ballpark  figure  for  determin¬ 
ing  the  Fried  Parameter  for  a  given  system.  Statistical  methods  may  not  be  the  most 
accurate  or  complete  descriptions  but  may  provide  a  good  enough  representation  for 
an  aero-optical  flow  over  an  Airborne  Laser  Weapons  System. 

Based  on  the  results  presented  in  Chapter  4,  this  particular  aero-optical  flow  (2-D 
shear  mixing  layer)  does  not  seem  to  have  Kolmogorov  like  defining  parameters  which 
could  be  a  result  of  not  capturing  the  full  distribution  of  eddy  scales.  The  structure 
functions  do  not  seem  to  behave  with  the  expected  trends.  The  structure  constant  is 
not  a  constant  and  it’s  range  of  values  is  higher  than  what  is  expected  for  the  turbulent 
atmosphere.  However  there  are  ranges  where  the  data  seems  to  reinforce  the  values 
derived.  This  gives  some  providence  in  investigating  deeper.  Some  opportunities  for 
improvement  this  are  found  in  the  next  section. 

5.2  Opportunities  for  Improvement 

The  first  opportunity  for  improvement  is  to  choose  one  index  of  refraction  calcu¬ 
lation.  The  results  show  that  there  is  very  little  difference  between  the  Edlen  and 
Gladstone-Dale  methods  at  the  near-IR  wavelengths.  In  light  of  this,  selecting  and 
using  Edlen  would  be  sufficient.  It  accounts  for  both  wavelength  and  density  varia¬ 
tions. 

The  second  opportunity  is  to  compare  multiple  reference  points  for  the  structure 
function  development.  This  ignores  the  feasibility  of  assuming  the  flow  is  isotropic 
but  may  provide  some  smoothing  to  the  structure  functions  that  may  make  them 
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more  Kolmogorov.  The  isotropic  nature  may  be  restored  if  these  various  reference 
point  parameters  are  averaged. 

The  third  opportunity  is  to  improve  the  resolution  of  the  mesh.  This  would  po¬ 
tentially  better  smooth  out  the  data  as  well  as  provide  more  data  points.  However, 
this  would  be  a  more  intensive  operation  computationally.  It  would  require  rebuilding 
the  fluid  grid  to  capture  more  flow  data,  which  would  increase  computational  require¬ 
ments  and  time.  The  additional  data  would  also  require  more  work  from  Matlab  to 
generate  the  meshes  from  which  the  parameters  are  developed  from. 

5.3  Future  Research  Opportunities 

This  investigation  was  confined  to  a  very  isolated  case.  A  2-D  shear  mixing  layer 
case  provides  a  good  opportunity  to  explore  the  details  of  the  methodology  developed 
for  this  research.  However,  it  is  apparent  that  there  are  opportunities  to  improve 
this  methodology.  Depending  on  the  feasibility  or  the  direction  of  implementing 
those  opportunities,  it  would  be  important  to  have  an  idea  of  how  to  extend  this 
methodology  for  more  rclevent  scenarios. 

The  next  logical  step  is  to  make  this  2-D  shear  mixing  layer  fluid  field  evolve  with 
time  and  what  effects  unsteadiness  in  the  flow  has  on  this  parameterization.  This 
can  begin  with  the  initial  solution  used  for  this  project  and  take  sample  data  from 
several  timesteps  as  it  evolves.  The  subsequent  logical  step  would  then  to  make  a  3-D 
laminar  flow  field.  The  more  interesting  case  would  then  be  a  3-D  turbulent  cases 
which  may  behave  more  inline  with  Kolmogorov  Theory.  The  problem  with  a  3-D 
turbulent  case  would  be  the  increase  in  computational  requirements.  The  2-D  grid 
for  this  research  included  500,000  grid  cells.  A  3-D  grid  would  have  n  x  500,  000  more 
cells  depending  on  the  minimum  gridspace  required.  That  number  would  likely  be 
much  higher  in  order  to  capture  the  boundary  layer  development  for  a  turbulent  case 
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as  this  grid  had  a  Y+  ten  times  greater  than  a  viable  Y+  for  a  turbulent  case. 

The  ultimate  goal  for  this  methodology  is  for  it  to  be  expanded  onto  a  turret 
configuration.  A  turret  would  provide  an  operational  configuration  for  developing 
aero-optical  turbulence  and  determining  the  Kolmogorov  relation  feasibility.  This 
would  also  introduce  a  lot  more  complexity  in  the  development  of  the  grid  and  com¬ 
putational  requirements  to  apply  this  methodology.  However,  it  would  make  a  more 
definitive  conclusion  regarding  the  Kolmogorov-like  relationship  of  an  aero-optical 
flow.  If  this  scenario  had  a  Kolmogorov-like  aero-optical  turbulence,  system  design 
could  be  simplified  by  applying  simple  statistical  methods  for  accounting  for  turbu¬ 
lence,  versus  extensive  simulations  of  chaotic  aero-optical  turbulence  fields. 
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Appendix  A.  Boundary  Condition  File 

################################################## 

Boundary  Condition  Specification  File  for: 

'Replace  this  line  with  a  title' 
################################################## 

9 

top_inlet 

Source 

Riemann  Invariant 

P-Stat  T-Stat  K  Omega  Mach  Axis  End  Points  Swirl  (A,B,C) 
101325.0  288.15  -1.  -1.  0.6  0.  0.  0.  1.  0.  0.  0.  0.  0. 

No 

################################################ 

10 
top 

Solid  Wall 

Slip 

No 

################################################ 

11 

outflow 
Farf ield 

Modified  Riemann  Invariant 

P-Stat  T-Stat  K  Omega  Mach  Alpha  Beta 

-1.  -1.  -1.  -1.  -1.  -370.  -370. 

################################################ 
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12 


bottom 
Solid  Wall 
Slip 
No 

################################################ 

13 

bottom_inlet 

Source 

Riemann  Invariant 

P-Stat  T-Stat  K  Omega  Mach  Axis  End  Points  Swirl  (A,B,C) 
101325.0  288.15  -1.  -1.  0.1  0.  0.  0.  1.  0.  0.  0.  0.  0. 

No 

################################################ 

14 

splitter_plate 
Solid  Wall 
Isothermal  No  Slip 
Wall  Temperature 
288.15 
No 

################################################ 
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Appendix  B.  AVUS  Job  File 


# ! /bin/bash 

######################### 

#  Queue  Options 
######################### 

#PBS  -1  nodes=6:ppn=2 

#PBS  -M  james.bowers@afit.edu 
#PBS  -m  abe 
#PBS  -j  oe 
#PBS  -V 

# 

#echo  Working  directory  is  $PFS_0_W0RKDIR 
cd  $PBS_0_W0RKDIR 
######################### 

#  Script  Banner 
######################### 
echo  -e  "" 

echo  -e  "======================" 

echo  -e  "  AVUS  Job  File  Script  " 
echo  -e  "======================" 

######################### 

#  AVUS  -  File  names 
######################### 

export  GRIDNAME=spltplt_2D_laminar ; 
export  BCNAME=spltplt_2D_laminar ; 
export  RESULTNAME=spltplt_2D_laminar ; 
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export  TAPNAME=; 

######################### 

#  AVUS  -  File  paths 
######################### 

#export  AVUSL0C=${H0ME}/ avus/bin ; 
export  AVUSL0C=${H0ME]-/ avus ; 

export  JOBLOC=${HOME}-/mydocs/runschtuf  f ; 
export  GRIDL0C=${J0BL0C}- ; 
export  BCL0C=${J0BL0C]- ; 

export  SCRATCH=${PBS_0_W0RKDIR}/ spltplt_2D_laminar ; 
export  RESULTLOC=${SCRATCH} ; 
export  TAPLOC=${SCRATCH} ; 

######################### 

#  AVUS  -  Executable  Spec 
######################### 

export  MACHINE_ARCH="linux" ;  #  linux  |  macosx  |  etc... 

export  PRECISIDN="double" ;  #  single  Idouble 

#export  RUNSCRIPT= " avus . linux. dp" ;  #  AvusRUN  |  AvusRUN_ibm 
export  RUNSCRIPT=" AvusRUN";  #  AvusRUN  |  AvusRUN.ibm 
######################### 

#  MPI  -  Run  command 
######################### 
export  RUN="mpirun" 

######################### 

#  Clean  scratch  directory 
cd  $SCRATCH/ 
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rra  -f  avus  AvusIN*  AvusOUT*  fort.*  *. shutdown 


######################### 

cat  >  $SCRATCH/$RESULTNAME. inp  «  EOF 

TITLE 

Shear  Flow  Turbulence  DES  24  Nov  2009 

INPUT  FILE  CONTROL  PARAMETERS 

START  OPTION  (1=INITIAL  RUN,  2=RESTART ,  3=RESTART  &  RECALC  WALL  DIST) 
1 

NO.  PROCESSORS  GRID  &  INTERSECTION  FILE  FORMAT  (0=UNF0RM,  1=F0RM) 
10  0 

SPLITTING  PROCS  PROVIDE  FLOW  DATA?  (0=N0 , 1-11=YES) 

-1  0 

OUTPUT  FILE  CONTROL  PARAMETERS 

CREATE  PICTURE  FILE?  (0=N0 , 1-9=YES)  FORMAT (0=UNF0RM,  1=F0RM) 


6 

CONVERGENCE  FREQ. 

0 

RESTART  FREQ. 

MOVIE  TAP  FREQ. 

PIX  FREQ 

5 

-1 

-1 

500 

ALGORITHM  PARAMETERS 
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EQUATION  SET  (1=EULER,  2=LAMINAR  N-S,  3=TURBULENT  N-S) 

2 

TURBULENCE  MODEL  (1=SPALART , 2=Spalart-DES , 3=MentorBSL , 4=MentorSSTT) 
2 

SPATIAL  ACCURACY  (1  OR  2)  TEMPORAL  ACCURACY  (1  OR  2) 

2  1 

THETA  (0. 5-1.0) 

1.00 

RHS  IFLUX(1=G&G)  LSTSQ . WTS (0=0FF , 1=0N)  LIMITER (0=0FF , 1=B&J , 2=Venk) 

1  0  2 

ITERATIVE  MATRIX  SOLUTION  SCHEME  (1=JAC0BI,  2=GAUSS-SEIDEL) 

2 

NO.  ITERATIONS  (SWEEPS)  OF  ABOVE  MATRIX  SOLUTION  SCHEME 
32 

INVISCID  JACOBIAN  DDF  VISCOUS  JACOBIAN  DDF 

0.1  0.1 

CFL  TIME  STEPS  NEWTON  SUB-ITERATIONS 

l.e+6  2000  1 

TIME  ACCURATE?  REQUESTED  TIME  STEP 

0  -1. 

REFERENCE  CONDITIONS  &  PHYSICAL  CONSTANTS 

UNITS  (1=MKS,  2=CGS ,  3=F00T-SLUG-SEC ,  4=INCH-SNAIL-SEC) 

1 

MACH  NO.  ANGLE  OF  ATTACK  ANGLE  OF  SIDESLIP 
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0. 


0. 


0. 


STATIC  PRESSURE  STATIC  TEMPERATURE 
101325.0  288.15 

GAMMA  GAS  CONSTANT  PRANDTL  NUMBER  GRAVITY 
-1.  -1.  -1.  0. 

INITIAL  CONDITIONS 

MACH  NO.  ANGLE  OF  ATTACK  ANGLE  OF  SIDESLIP 
0.  0.  0. 

STATIC  PRESSURE  STATIC  TEMPERATURE 
101325.0  288.15 

TURBULENT  KINEMATIC  VISCOSITY  RATIO 
-1.  -1. 

GEOMETRY  PARAMETERS 

COORDINATE  SYSTEM  (1=FL057,  2=PANAIR,  3=AXI -SYMMETRIC) 
1 

AXISYMMETRIC  FORCE  ACCOUNTING 
-1 

REFERENCE  AREA 

0.0 

X,Y,Z  COORDINATES  OF  MOMENT  REFERENCE  POINT 
0.0  0.0  0.0 

REFERENCE  LENGTHS  FOR  MOMENTS  ABOUT  X-,Y-  AND  Z-AXIS 
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0.0  0.0  0.0 


NONINERTIAL  REFERENCE  FRAME 

OMEGA (1 : 3)  -  NON-INERTIAL  ROTATION  VECTOR  (RAD/S) 

0.00  0.00  0.00 

0MEGACNTR(1 : 3)  -  CENTER  OF  NON-INERTIAL  ROTATION 
0.00  0.00  0.00 

INITIAL  ROTATION  (0=N0NE, 1=USE  OMEGA, 2=USE  INITOMEGA) 

0 

INITOMEGA (1 : 3)  -  INITIAL  NON-INERTIAL  ROTATION  (RAD/S) 
0.00  0.00  0.00 

INIT0MEGACNTR(1 : 3)  -  CENTER  OF  INITIAL  ROTATION 
0.0  0.0  0.0 

FLUID  VOLUME  MESH  DEFORMATION 

MESH  DEFORMATION  METHOD  (0=N0NE,  1=V0LDEF,  2=DYNMESH) 

0 

END  OF  INPUT  INFORMATION 

EOF 

# - 

#  AvusRUN  script  argument  list: 
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# 

#  1  -  AVUS  input  file  name 

#  2  -  AVUS  output  file  name 

#  3  -  precision  switch 

#  4  -  grid  file  name 

#  5  -  old  restart  file  name 

#  6  -  new  restart  file  name 

#  7  -  picture  file  name 

#  8  -  tap  location  file  name 

#  9  -  shutdown  file  name 

#  10  -  performance  file 

#  11  -  bl  trip  file  name 

#  12  -  be  file  name 

#  13  -  overwrite  flag 

#  14  -  machine  type 

#  15  -  J -mdiceargs’  (optional : only  when  run  from  MDICE) 

#  16  -  string  of  mdice  args  (optional : only  when  run  from  MDICE) 

#  - 

#  Available  precision: 

#  single,  double 

# 

#  Available  machine  types : 

#  ibm,sgi,t3e 

#  - 

# 

# 
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$AVUSLOC/$RUNSCRIPT  \ 
$SCRATCH/$RESULTNAME. inp  \ 
$RESULTLOC/$RESULTNAME. out  \ 
$PRECISION  \ 

$GRIDL0C/ $GRIDNAME . grd  \ 
$SCRATCH/junk. intr  \ 
$RESULTLOC/$RESULTNAME.rst  \ 
$SCRATCH/ $RESULTNAME . tr st  \ 
$RESULTLOC/$RESULTNAME.pix  \ 
$SCRATCH/$ JOBNAME . tap  \ 
$SCRATCH/$RESULTNAME. shutdown  \ 
$SCRATCH/$RESULTNAME . movtap  \ 
$SCRATCH/ramp . trip  \ 
$BCLOC/$BCNAME . be  \ 
overwrite  \ 

$MACHINE_ARCH  \ 
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Appendix  C.  Parameterization  Algorithm 


Listing  C.l.  flowoptics.m 


Wo  FlowOptics 

%Written  by:  Captain  James  C  Bowers 
%Date:  9  Dec  2009 

^Summary :  This  code  uses  flow  data  generated  using  a  CFD  solver  and 

%calculates  Statistical  Optical  Parameters 


Wo  Data  Import 
dat  areadin 

Wo  Pressure  Temperature  Conversion 
ptconv 

Wo  Index  of  Refraction  Calculation 
index 

Wo  Allocating  data  to  Grids 
datamesh 

Wo  Radius  Calculation  Mesh 
radmesh 

Wo  OPL  Calculation  Loop 
oplcalc 

Wo  OPD  Calculation  Loop 
opdcalc 

Wo  Theta  &  Amplitude  calc 
phaseamp 

Wo  Developing  Structure  Functions 
structure 

Wo  Flowfield  Data  Plotting  for  Visualization 
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flow  viz 


%Zo  Data  Presentation 
metric 


%  End  of  Routine 
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Listing  C.2.  datareadin.m 


% - dataread.m - % 

%Reads  in  result  data  output  by  blacksmith  for  given  flow  data 
%Created  by:  James  C  Bowers  Date:  5  Jan  2010 

%Used  for  pulling  the  flowdata  into  arrays  for  calculations  by  itself  or  as 
%part  of  a  subroutine  in  larger  code  to  analyze  data  for  optical 
%properties.  jcb  5  jan  2010 

% - % 

%Zo  Clears  and  Closes  any  open  files  and  variable  used  or  created  in  routine 
close  all 

clear  fid  i  tline  numptsstr  numpts  x  y  rho  u  v  w  press  kmu 

Wo  Opens  Data  File  for  read  in  and  metrics 

f  id=fopen  (  1  /  home/  af  i  t  en  3  /  gaplOm  /jcbowers  /  Mat  lab  /data/2D_laminar_data.dat  '  ,  'r  ' 
i  =0; 

%Zo  Number  of  Points  of  Data  to  input 

%While  loop  goes  line  by  line  through  header  to  find  the  spot  in  the  third 
%line  of  the  header  which  tells  the  number  of  data  points  to  pull  in  since 
%there  is  more  lines  of  data  regarding  the  geometry  past  the  area  of  data 
%that  is  of  interest. 

while  feoff  fid)==0 

tline=fgetl  (fid  ); 
i=i  +1; 

%When  it  hits  the  third  line  it  pulls  the  string  with  the  number  of 
%data  points  and  then  converts  it  to  a  number.  This  method  is  specific 
%to  the  file  I  am  reading  in.  May  need  to  be  fixed  or  adjusted  with 
%different  files. 

i  f  i  ==3 

numptsstr=t line  (26:31); 
numpts=st  r2double  (  numptsstr  )  ; 
break 


end 
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end 


%Closes  file  after  getting  numpts  for  read 
fclose  (  fid  )  ; 

%/o  Data  Read  In 

%Ignores  first  3  lines  of  data  file  then  stores  the  following  tabular  data 
%into  respective  arrays  for  the  number  of  points  (numpts)  it  should  get. 

[x  ,  y  ,  rho  ,u,v,w,  press  ,  kmu]  =  text  read  .  .  . 

(  1  /  home/  afi  t  en3  /  gap  10m/ j  cbowers  /  Mat  lab  /data/2D_laminar_data.dat  1  .  .  . 

,  '%f  %f  %f  %f  %f  %f  %f  %f  \n  '  , numpts,  'headerlines  '  ,3); 

Wo  Data  Clean  Up 

%Not  all  of  the  flow  data  is  useful  and  will  be  cleared  to  free  of  space 
%and  clutter  of  workspace. 

clear  u  v  w  kmu  tline  i  fid  ans  numptsstr 

%  End  of  Routine 
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Listing  C.3.  ptconv.m 


% - ptconv.m - 

%Reads  in  pressure  and  density  data  to  create  temperature  data  and  then 
Reconverts  the  pressure  from  Pascals  to  mbars. 

%Created  by  James  C  Bowers  Date:  5  Jan  2010 
% - 

%Inputs  : 

%numpts  =  numpts  of  data  from  file 
%press  =  pressure  data  from  file 
%Outputs  : 

%temp  =  temperature  data  calculated 

% - 

clear  temp  i 

temp=zeros  ( numpts  ,  1 )  ;  %I  n  i  t  i  a  1  i  z  i  n  g  temperature  array 

for  i  =1:  numpts 

temp  ( i  ,  1 )  =  press  ( i  ,  1 )  /  (  rho  ( i  ,  1 )  *  2  8  7 ) ;  %Temperature  calculation 
press  ( i  ,1)  =  press  ( i  ,  1)  /  1 0 0  ;  RoPressure  Conversion 


end 


clear  i 


%  End  Routine 


Listing  C.4.  index. m 


% - index,  m - % 

%Routine  that  take  flow  data  and  calculates  index  of  refraction  from  the 
%data  at  those  points. 

%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%Input  Variables 

%press  =  pressure  data  (called  in) 

%temp  =  temperature  data  (called  in) 

%rho  =  density  data  (called  in) 

%numpts  =  number  of  data  points  (called  in) 

%lamda  =  wavelength  in  microns  (first  specified) 

%kgd  =  Gladstone— Dale  coefficient  (first  specified) 

% 

%Output  Variables 

%ned  =  index  of  refraction  using  Edlen 

%ngd  =  index  of  refraction  using  Gladstone— Dale  (  density ) 

% - % 

clear  ned  ngd  lamda  kgd  i 

%I n i t i al i z i n g  refractive  index  grids 
ned=zeros  (numpts  ,  1 )  ; 
ngd=zeros  (numpts  ,  1 )  ; 

lamda  =  l;  %  Wavelength  in  microns 

kgd=2.25e—  4;  %Gladstone  Dale  Constant  assuming  1  micron  wavelength 

for  i  =1:  numpts 

%Edlen  Index  of  Refraction 

ned  ( i  ,  1  )  =  1  +  ( 10  ~  (- 8)  )*  (8342 . 1 2  +2406030/(130  - 1/ (lamda  ~  2 ) )  .  .  . 

+  (15997/(38.9  —  1/lamda  ~2)))*(rho(i  ,  1 )  /  1  .221  )  ; 


%Calculates  index  of  refraction  using  Gladstone— Dale  for  1  micron 
%wavelength  and  the  density  of  the  flowfield 
ngd  ( i  ,l)  =  l  +  kgd*rho  ( i  ,1); 
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end 


%Clean  Up 
clear  kgd  i 

%end  of  Routine 
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Listing  C.5.  datamesh.m 


% - datamesh.m - % 

%Routine  that  takes  flow  data  and  puts  the  data  into  a  grid  for  iterating 
%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%Input  Variables 
%x=  xdata  (called  in) 

%y=  y  data  (called  in) 

%press  =  pressure  data  (called  in) 

%temp  =  temperature  data  (called  in) 

%rho  =  density  data  (called  in) 

%ned  =  index  of  refraction  using  edlen  (called  in) 

%ngd  =  index  of  refraction  using  Gladstone— Dale  ( density )  (called  in) 

%xx  =  number  of  grid  points  in  x  direction  (first  specified) 

%yy  =  number  of  grid  points  in  y  direction  (first  specified) 

% 

%Output  Variables 

%xmid  =  grid  midpoint  in  x  direction 

%ymid  =  grid  midpoint  in  y  direction 

%xmat=  x  point  array 

%ymat  =  y  point  array 

%tempgrd  =  Temperature  Grid 

%pressgrd  =  Pressure  Grid 

%rhogrd  =  Density  Grid 

%xgrd  =  X  Grid 

%ygrd  =  Y  Grid 

%nedgrd  =  Grids  the  refractive  index  from  edlen 
%nedgrd  =  Grids  the  refractive  index  from  G— D 
% - % 

%Pre— Clean  Up 

clear  xx  yy  xmid  ymid  xmat  ymat  tempgrd  pressgrd  rhogrd  xgrd  ygrd  nptgrd  ngdgrd 

%Dimension  Setup 
xx  =  1001 ; 
yy — 1001 ; 

%Grid  Midpoints 
xmid=c  e  i  1  (xx  /2); 
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ymid=c  e  i  1  (yy  /  2) ; 


%Below  is  Location  allocation  for  gridding  the  data.  This  will  eventually 
%change  with  the  next  iteration  be  conditional  upon  either  data  a/o  path 
xmat=linspace  (0 , 1  ,xx  ) ;  %X  Location  Allocation 
ymat=linspace  (0 .5  ,  —  0 .5  ,  yy  )  Location  Allocation 

%Setting  up  Grids  for  current  data 

tempgrd=gr  iddat  a  (x  ,  y  ,  temp  ,  xmat  ,  ymat )  ;  %Temperature  Grid 
pressgrd=gr iddat  a  (x  ,  y  ,  press  ,  xmat ,  ymat ) ;  %Pressure  Grid 
rhogrd=gr  iddat  a  (x  ,  y  ,  rho  ,  xmat ,  ymat )  ;  %Density  Grid 
xgrd=griddat  a  (x  ,  y  ,  x  ,  xmat ,  ymat )  ;  %K  Grid 
ygrd=griddat  a  (x  ,  y  ,  y  ,  xmat ,  ymat )  ;  Grid 

nedgrd=gr  iddat  a  (x  ,  y  ,  ned  ,  xmat ,  ymat )  ;  %Grids  the  refractive  index  from  Edlen 
ngdgrd=gr  iddat  a  (x  ,  y  ,  ngd  ,  xmat ,  ymat )  ;  %Grids  the  refractive  index  from  G— D 

%End  of  Routine 
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Listing  C.6.  radmesh.m 


% - radmesh.m - % 

%Routine  that  creates  radius  spread  from  x  &  y  points  in  gridded  format 
%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%Input  Variables 

%xgrd  =  X  Grid  (called  in) 

%ygrd  =  Y  Grid  (called  in) 

%Output  Variables 

%ngrd  =  Grids  the  radius  from  x  &;  y  grids 
% - % 

clear  rgrd  i  j 

rgrd=zeros  (xx  ,  yy  )  ;  Reinitializing  radius  array 

for  i  =  l:xx 

for  j=l:yy 


%cleans  up  the  left  hand  column  of  x  points 
if  i==l 

xgrd(j  ,  1 )  =  0 ; 

end 

%Calculates  radius  from  central  point  in  flow  field  for  structure 
%function  calculations 

rgrd  ( i  ,  j)=sqrt  ( ( (  xgrd  ( i  ,  j  )-xgrd  (xmid  ,  ymid ) )"  2 )  +  ( ( ygrd  ( i  ,  j)-ygrd  (xmid  ,  ymid  ) )  *  2) ) ; 


end 

end 

clear  i  j 

%end  of  routine 
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Listing  C.7.  oplcalc.m 


% - o  pic  ale .  m - % 

%Routine  that  calculates  optical  path  length 
%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%Input  Variables 

%nedgrd  =  Grids  the  refractive  index  from  edlen  (called  in) 

%ngdgrd  =  Grids  the  refractive  index  from  G— D  (called  in) 

%dy  =  y  gridspacing  (first  specified) 

%suminted  =  integration  sum  for  edlen 
%sumintgd  =  integration  sum  for  gladstone  — dale 
%xx  ,  yy  =  grid  dimensions 
%Output  variables 

%oplgdmean  =  optical  path  length  mean  for  row  using  gladstone —dale 
%opledmean  =  optical  path  length  mean  for  row  using  edlen 
%opled  =  optical  path  length  integrated  up  the  column  using  edlen 
%oplgd  =  optical  path  length  integrated  up  the  column  using  glad  — dale 
% - % 

%Pre— Clean  up 

clear  dy  suminted  sumintgd  opled  oplgd  oplgdmean  opledmean  i  j 

dy=0.001;  %Delta  y  —  used  for  intergration  —  comes  from  length  and  spacing 

%Initializing  integration  sum 
suminted  =0; 
sumintgd  =0; 

%i  n  i  t  i  a  1  i  z  i  n  g  Matricies 
opled=zeros  (yy  ,  xx  )  ; 
oplgd=zeros (yy , xx ) ; 
oplgdmean=zeros  (yy  ,  1 )  ; 
opledmean=zeros  (yy  ,  1 )  ; 

for  j  =1  :xx  %Steps  through  columns  for  x  direction 

for  i=yy:— 1:1  %Starts  at  bottom  row  of  column  and  works  up 
%Integrates  up  column 
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suminted=suminted+nedgrd  ( i  ,  j  )*dy  ; 
sumintgd=sumintgd+ngdgrd  ( i  ,  j  )*dy  ; 

%S tores  OPL  at  those  points 
opled(i  ,  j  )  =  suminted  ; 
oplgd  ( i  ,  j  )  =  sumintgd  ; 


end 

%Re  sets  sumints 
suminted  =0; 
sumintgd  =  0; 


end 

%Calculates  the  mean  of  the  OPL' s  as  it  rises  from  the  bottom  by  row 
for  i=yy :  —  1 : 1 

opledmean  ( i  ,  1  )  =  mean  (opled(i  ,  :  )  )  ; 
oplgdmean  ( i  ,  1  )  =  mean  (oplgd(i 

end 

%Post  Clean  Up 

clear  suminted  sumintgd  i  j 
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Listing  C.8.  opdcalc.m 


% - opdcalc.m - % 

%Routine  that  calculates  optical  path  difference 
%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%Input  variables 

%xx  ,  yy  =  grid  dimensions 

%oplgdmean  =  optical  path  length  mean  for  row  using  gladstone  — dale 
%opledmean  =  optical  path  length  mean  for  row  using  edlen 
%opled  =  optical  path  length  integrated  up  the  column  using  edlen 

%oplgd  =  optical  path  length  integrated  up  the  column  using  G— D 

%Output  Variables 

%opded  =  optical  path  difference  using  press— temp 

%opdgd  =  optical  path  difference  using  press— temp 

% - % 

clear  i  j  opded  opdgd 

opded=zeros (yy , xx ) ; 
opdgd=zeros (yy , xx ) ; 

for  j  =l:xx 

for  i=yy :  —  1 : 1 

opded  ( i  ,  j  )  =  opled  ( i  ,  j )  — opledmean  ( i  ,  1 )  ; 
opdgd  ( i  ,  j  )=oplgd  ( i  ,  j  )-oplgdmean  ( i  ,  1 )  ; 


end 


end 


clear  i  j 
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Listing  C.9.  phaseamp.m 


% - phaseamp.m - % 

%Routine  that  calculates  the  phase  from  opd  and  amplitude  of  the  beam  for 
%use  in  subsequent  metric  calculations. 

%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%Input  variables 

%xx  ,  yy  =  grid  dimensions 

%opdgd  =  optical  path  difference  using  gladst  one —dale  (called— in) 

%opded  =  optical  path  difference  using  edlen  (called— in) 

%wo  =  beam  waist  radius  provided  (first  specified) 

%lamda  =  beam  wavelength  (called  in) 

%Output  Variables 

%thetaed  =  beam  phase  at  points  using  edlen 

%thetagd  =  beam  phase  at  points  using  G— D 

%a  =  beam  amplitude  at  points  (assumes  normalized  peak  of  1  at  center) 

% - % 

clear  wo  thetapt  thetagd  i  j  a 

%Specifies  a  beam  waist  radius 
wo=0 . 5  ; 

%l  n  i  t  i  al  i  z  i  n  g  Output  matricies  for  speed 
thetaed=zeros  (xx  ,  yy  ) ; 
thetagd=zeros  (xx  ,  yy  ) ; 
a=zeros (xx , yy  ) ; 

%Phase  and  Amplitude  Calculation  Iteration 
for  i=l:xx 

for  j=l:yy 


%Phase  at  points  using  OPD  Sz  wavelength 
thetaed(i  ,  j  )  =  ((—  2*pi)/(  lamda*10~  —  6))*opded  ( i  ,  j  ) ; 
thetagdfi  ,  j  )  =  ((—  2*  pi )  /  ( lamda*10/'  —  6))*opdgd  ( i  ,  j  ); 

%Natural  Log  of  amplitude  assuming  Gaussian  beam 
a(i  ,  j  )  =  -((xgrd  ( i  ,  j  )-xgrd  (xmid  ,ymid) )  ‘  2)  /  (wo"  2) ; 
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end 


end 


%Post  Clean  Up 
clear  i  j 


%End  of  Routine 
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Listing  C.10.  structure. m 


%- 


-structure.m- 


-% 


%Routine  that  calculates  the  structure  functions  and  then  determinens  the 

%Cn''2's  from  various  methods 

%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%/o  Input  Variables 

%dr  =  radial  step  size  (first  specified) 

%tol=  radius  incorporation  tolerance  for  data  points  (first  specified) 
%nedgrd  =  Grids  the  refractive  index  from  Edlen  (called  in) 

%ngdgrd  =  Grids  the  refractive  index  from  G— D  (called  in) 

%xx  ,  yy  =  grid  dimensions 
%xmid  ,  ymid  =  grid  midpoints 

%rgrd  =  flow  radius  from  center  point  (called  in) 

%tempgrd  =  temperature  data  from  flow  (called  in) 

%thetapt  =  beam  phase  at  points  using  press— temp  (called  in) 

%thetagd  =  beam  phase  at  points  using  press— temp  (called  in) 

%a  =  beam  amplitude  at  points  (assumes  normalized  peak  of  1  at  center) 
Output  variables 

%nmax  =  calculated  max  number  of  radial  steps 


of 

index  difference  for 

average  from 

G-D 

of 

index  difference  for 

average  from 

Edlen 

sum 

of  phase  difference 

for  average 

from 

Edlen 

sum 

of  phase  difference 

for  average 

from 

G-D 

%dasum  =  sum  of  amplitude  difference  for  average 

%dtsum  =  sum  of  temperature  difference  for  average 

%k  =  number  of  data  points  at  radius  used  for  average 

%pntcnt  =  collected  number  of  data  points  for  all  radii 

%r  =  radius  values 

%dned  =  Structure  function  edlen 

%dngd  =  Structure  Function  G— D 

%cn2ed  =  Cn2  calculated  from  edlen 

%cn2gd  =  Cn2  calculated  from  G— D 

%dthetaed  =  Structure  Function  for  phase  for  Edlen 
%dthetagd  =  Structure  Function  for  phase  for  G— D 
%da  =  Structure  Function  for  In  ( amplitude ) 

%dtaed  =  Combined  Struct  Function  for  Amp  &  Phase  for  Edlen 
%dtagd  =  Combined  Struct  Function  for  Amp  &  Phase  for  G— D 
%dt  =  Struct  Function  for  Temperature 
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%ct2  =  Ct2 

%cn2t  =  Cn2  calculated  from  Ct2 
%roed  =  fried  coherance  length  from  Edlen 
%rogd  =  fried  coherance  length  from  G— D 
% - 


Pre  Clean 

clear  nmax  dr  tol  nmax  dngdsum  dnedsum  dasum  dtsum  k  i  j  n 
clear  pntcnt  r  dned  dngd  cn2ed  cn2gd  dthatapt  dtheatgd 
clear  da  dtaed  dtagd  dt  ct2  cn2t  roed  rogd 

nmax=floor  (xx/2) ;  %Limits  the  number  of  radial  steps  out  from  the  center 
dr=0.001;  %Radial  step  size 

%  Because  grid  is  cartesian  and  the  calculations  are  based  on  polar 
%  coordinates  I  specify  a  radius  tolerance  to  smooth  out  the  number  of 
%  sample  points  within  my  radial  grid 
tol  =0 .0005  ; 


-% 


%  Initializing  sums  and  counters  prior  to  iteration 

dnedsum=0; 

dngdsum=0; 

dthetaedsum  =  0; 

dthetagdsum  =  0; 

dasum  =  0; 

dtsum  =  0; 

k=0;  %Sample  point  counter 
r=zeros  ( 1  ,nmax) ; 
dned=zeros  (1  ,nmax)  ; 
dngd=zeros  (1  ,nmax)  ; 
cn2ed=zeros  ( 1  , nmax )  ; 
cn2gd=zeros  (1  ,nmax) ; 
dthet aed=zer os  ( 1  , nmax ) ; 
dthetagd=zer os  ( 1  , nmax ) ; 
da=zeros  ( 1  , nmax )  ; 
dtaed=zeros  (1  ,nmax) ; 
dtagd=zeros  ( 1  , nmax )  ; 
dt=zeros  ( 1  , nmax )  ; 
ct2=zeros  ( 1  ,nmax )  ; 
cn2t=zeros  (1  ,nmax)  ; 
roed=zeros  (1  ,nmax)  ; 
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rogd=zeros  (1  ,nmax)  ; 
ptcnt=zeros  (1  ,nmax) ; 

for  n  =  l:nmax  %Radial  outward  step 

r(n)=n*dr;  %Stores  current  radius 

for  i  =  1 :  xx 

for  j  =  1 : y y 


%Steps  through  the  radius  grid  to  see  if  the  radius  point  data 
%is  with  in  the  step  size  specified  radius  +/—  the  tolerance. 

%If  radius  meets  those  criteria  it  uses  the  data  point  to 
%calculate  the  structure  functions 
if  rgrd  ( i  ,  j  )<r  (n)+ 1 ol  &&;  rgrd ( i  , j )>r (n)— t ol 

%Refractive  Index  Structure  Function 

dnedsum=dnedsum+(abs  ( (  nedgrd  ( i  ,  j )  — nedgrd  (xmid  ,  ymid  ) ) )  "  2  ) ; 

dngdsum=dngdsum+(abs  ( (  ngdgrd  ( i  ,  j )  — ngdgrd  (xmid  ,  ymid  ) ) )  ~  2  ) ; 

%Temperature  Structure  Function 

dtsum=dtsum+(abs  ( ( tempgrd  ( i  ,  j )  —  tempgrd  (xmid  ,  ymid  ) ) )  "  2  ) ; 

%Phase  Structure  Function 

dthetaedsum=dthetaedsum+(abs  ((thetaed(i  ,  j)  —  thetaed  (xmid  ,  ymid  ) )  "  2  ) ) ; 

dthetagdsum=dthetagdsum  +  (abs  ((thetagd(i  ,  j)  —  thetagd  (xmid  ,  ymid  ) )  "  2  ) ) ; 

/^Amplitude  Structure  Function 

dasum=dasum+abs  ( (  a  ( i  ,  j )  — a  (xmid  ,  ymid  ) )  "  2  ) ; 

k=k  +  l;  %Sample  Count  Number  for  averaging 


end 


end 

end 


%Average  of  Dn  for  the  specified  radius 
dned  (n)=dnedsum/k ; 
dngd  (n)=dngdsum/k ; 
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%Calculates  Cn"2  from  Refractive  Index  Struture  function  assuming  rA(2/3) 
cn2ed  (n)  =  (dned  (n)  /((r(n))A(2/3))); 
cn2gd(n)=(dngd(n)/((r (n)) " (2/3))); 

%Calculates  Cn~2  from  Temperature  Struture  function  assuming  rA(2/3) 
dt  (n)=dtsum/k ; 

ct2  (n)  =  (dt(n)  /  ( (  r  (n ) )  ‘  (2/3))); 

cn2t  (n)  =  (ct2(n))*((79*10A( -6))*1001  .3/(288.15  *2)) ‘2; 

%Structure  function  for  Phase 
dthetaed  (n)  =  dthetaedsum /k ; 
dthetagd  (n)  =  dthetagdsum /k ; 

%Structure  function  for  Amplitude 
da(n)=dasum/k ; 

%Combined  Structure  Function 
dtaed  (n)=da (n)  + dthetaed  (n  )  ; 
dtagd (n)=da(n)+ dthetagd (n ) ; 

%Fried  Coherance  length  calculation 
roed(n)=r(n)/ (( dtaed (n)/6.88)A(3/5)); 
rogd (n)=r(n)/(( dtagd (n)/6 .88 ) ' (3/5)); 

%Calculates  Cn"2  from  fried  parameter 

cn2roed(n)  =  ((0.185  )*(((  lamda  *  1 0  A  (  —  6))A(6/5))/  ((lA(3/5))*(roed(n)))))  A  (5/3) 
cn2rogd(n)  =  ((0.185  )*(((  lamda  *  1 0  A  (  —  6))A(6/5))/  ((lA(3/5))*(rogd(n)))))  A  (5/3) 

%Calculates  ro  from  Cn~2 

rocn2ed(n)  =  (0.185  )*(((  lamda  *  1 0  A  (  —  6))  ^  (6/5))  /  ((lA(3/5))*(cn2ed(n)  A  (3/5)))); 
rocn2gd(n)=(0.185  )*(((  lamda  *  1 0  A  (  —  6))A(6/5))/  ((lA(3/5))*(cn2gd(n)  A  (3/5)))); 

%Calculates  Phase  Structure  function  from  Cn~2 

dthetacn2ed  (n)=2 . 9  1  *(((2*  pi )  /  ( lamda  *  1 0  A  (  —  6)))/'(2))*(r(n)~(5/3))*cn2ed(n); 
dthetacn2gd  (n)=2 . 9 1  *(((2*  pi )  /  ( lamda  *  1 0  A  (  —  6)))~(2))*(r(n)A(5/3))*  cn2gd  (n  ) ; 

%Stores  the  number  of  sample  points  for  pure  statistical  info 
ptcnt  (n)=k ; 
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%Resets  all  counters /sums  for  the  next  iteration 

dnedsum  =  0; 

dngdsum  =  0; 

dthetaedsum  =  0; 

dthetagdsum  =  0; 

dasum  =  0; 

k  =  0; 


end 

%/o  Post  Clean 

clear  i  j  k  dasum  dthetagdsum  dthetaedsum  dngdsum  dnedsum 
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Listing  C.ll.  flowviz.m 


% 

-% 


% 

clear  figcnt 

figcnt  =1; 

figure  (  figcnt  ) 

subplot  (1  ,3 ,1)  ,  mesh  ( xgrd  ,  ygrd  ,  tempgrd ) 
axis  ([0  1  -0.5  0.5  ]) 

t  i  1 1  e  (  '  Temperture  of  Flow  Field  (K)  '  ) 
xlabel (  'X'  ) 
y label (  'Y'  ) 

z label  (  '  Temperature  (K)  1  ) 

subplot  (1  ,3  ,2)  ,  mesh ( xgrd  ,  ygrd  ,  rhogrd  ) 

axis  ([0  1  -0.5  0.5  ]) 

t  i  1 1  e  (' Density  of  Flow  Field  (kg/nY3)') 
xlabel (  'X 1  ) 
y label (  ’Y  ' ) 

z  label  (' Density  (kg/m~3)') 

subplot  (1  ,3  ,3)  ,  mesh ( xgrd  ,  ygrd  ,  pressgrd) 

axis  ([0  1  -0.5  0.5  ]) 

t  i  1 1  e  (' Pressure  of  Flow  Field  (mbars)') 
xlabel (  'X 1  ) 
y label (  ’Y ' ) 

zlabel(  'Pressure  (mbars)  ') 


% - flow  viz .  m - 

%Routine  that  takes  flow  data  and  plots  it  for  visualization 
%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - 

%Input  Variables 

%tempgrd  =  Temperature  Grid  (called  in) 

%pressgrd  =  Pressure  Grid  (called  in) 

%rhogrd  =  Density  Grid  (called  in) 

%xgrd  =  X  Grid  (called  in) 

%ygrd  =  Y  Grid  (called  in) 

%nedgrd  =  Grids  the  refractive  index  from  Edlen  (called  in) 
%ngdgrd  =  Grids  the  refractive  index  from  G— D  (called  in) 
%figcnt  =  figure  counter  (first  specified) 

% - 
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figcnt  =  figcnt+l; 


figure  (  figcnt  ) 

subplot  (2 ,2  ,1)  ,  mesh ( xgrd  ,  ygrd  ,  nedgrd  ) 
axis  ([0  1  -0.5  0.5  ]) 

t  it  le  (  1  Refractive  Index(Edlen)  of  Flow  Field') 
xlabel( 'X') 
y label (  ' Y '  ) 
zlabel (  'n  '  ) 

subplot  (2 ,2  ,2)  ,  mesh  ( xgrd  ,  ygrd  ,  ngdgrd  ) 
axis  ([0  1  -0.5  0.5  ]  ) 

t  i  1 1  e  (  '  Refr  ac  t  i  ve  Index  (  Gladstone  Dale)  of  Flow  Field') 
xlabel( 'X') 
y label (  ' Y '  ) 
z labe 1 (  'n') 

subplot  (2 ,2  ,3)  ,  mesh  ( xgrd  ,  ygrd  ,  opdgd ) 
axis  ([0  1  -0.5  0.5  ]  ) 

t  i  1 1  e  (' Optical  Path  Difference  (Gladstone  Dale)  of  Flow  Field') 
xlabel (  'X'  ) 
y label (  ' Y '  ) 
zlabel (  'OPD  '  ) 

subplot  (2 ,2  ,4)  ,  mesh  ( xgrd  ,  ygrd  ,  opded  ) 
axis  ([0  1  -0.5  0.5  ]) 

t  i  1 1  e  (' Optical  Path  Difference  (Edlen)  of  Flow  Field') 
xlabel (  'X'  ) 
y label (  ' Y '  ) 
zlabel (  'OPD' ) 

figcnt  =  figcnt +1;  %updates  figure  counter 

%end  of  subroutine 
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Listing  C.12.  metric. m 


% - metric.m - % 

%Routine  that  takes  the  import  optics  data  calculated  and  plots  or  presents 

%that  data  in  a  visual  format 

%Created  by  James  C  Bowers  Date:  5  Jan  2010 

% - variables - % 

%Zo  Input  Variables 

%figcnt  =  figure  counter  (called  in) 

%opded  =  optical  path  difference  using  edlen  (called  in) 

%opdgd  =  optical  path  difference  using  G— D  (called  in) 

%xmat  =  linear  x  location  matrix  (called  in) 

%r  =  radii  for  basis  of  structure  funtion  from  center  (called  in) 

%dned  =  Structure  function  from  edlen  (called  in) 

%dngd  =  Structure  function  from  G— D  (called  in) 

%cn2gd  =  Cn~ 2  from  G— D  (called  in) 

%cn2edlen  =  Cn~2  from  edlen  (called  in) 

%cn2t  =  Cn~2  from  Temperature  Cn~2  (called  in) 

%roed  =  fried  coherence  length  using  edlen  (called  in) 

%rogd  =  fried  coherence  length  using  G— D  (called  in) 

Output  Variables 

%diff  =  absolute  difference  between  opd  '  s 
% - % 


Comparison  Plot 

%Plots  the  OPD' s  for  the  top  row  of  the  gridpoints 
d i f  f (1 , : )  =  abs ( opded ( 1 , : )  —  opdgd  (1  , : ) )  ; 
figure  (  figcnt  ) 

subplot  (2  ,1  ,1)  ,  plot  (xmat ,  opded  (1  , : )  , xmat ,  opdgd  (1  , : ) ) 

t  i  1 1  e  (  '  Opt  ical  Path  Difference  Along  Streamwise  Direction  Comparison') 
xlabel  (  'X(m)  '  ) 
y label  (  'OPD'  ) 

legend(  'Edlen'  ,  '  Gladstone— Dale  '  ,  'location  '  ,  'bestoutside  ') 
subplot  (2  ,1  ,2)  ,  plot  ( xmat ,  d i f f  ,  '  r  '  ) 

t  i  1 1  e  (' Opt  ical  Path  Difference  Comparison  Absolute  Difference') 
xlabel  (  'X(m)  '  ) 

y  label  (' Absolute  Difference') 


figcnt  =  figcnt+l; 


Data  Visualization 

figure  (  figcnt  ) 

loglog (r , dned , 'o' , r , dngd , ' xk 1 ) 
t  i  1 1  e  ('  Calculated  Structure  Function  with  r') 
x label (  ' r  (m)  '  ) 
y label (  'Dn'  ) 

legend(  'Edlen'  ,  1  Gladstone— Dale  '  ,  'Location  '  ,  'Bestoutside  ') 

figcnt=f  igcnt  +1; 

figure  (  figcnt  ) 

loglog  (r  ,  cn2ed  ,  1  og  '  ,  r  ,  cn2gd  ,  1  xk  '  ,  r  ,  cn2t  ,  '+r  '  ,r  ,cn2roed  ,  '  *m'  ,  r  ,  cn2rogd  ,  '  . b  '  ) 
t  i  1 1  e  (' Calculated  Cn~2  with  r 2  ~  ~  3  '  ) 
x label (  ' r  (m)  '  ) 
y 1 a  b  e  1  (  ' Cn "  2  '  ) 

legend(  'Edlen'  ,  'G— D '  ,  'CT"2  '  ,  '  r  _o  ( Edlen)  '  ,  '  r_o  (G— D)  '  ,  'Location  '  ,  'Bestoutside  ') 

figcnt=f  igcnt  +1; 

figure  (  figcnt  ) 

loglog (r  ,  roed  ,  'o'  ,  r , rogd  ,  ' xk  '  ,r  ,rocn2ed  ,  '  *m'  ,  r  ,  rocn2gd  ,  '  . b  '  ) 
t  i  1 1  e  (' Calculated  r_o') 
x label (  ' r  (m)  '  ) 
y label (  ' r _o  (m)  '  ) 

legend(  'Edlen'  ,  '  Gladstone— Dale  '  ,  'Cn~2  (Edlen )  '  ,  'Cn"2  (G— D)  '  ,  'Location  '  ,  'Bestoutside  ') 

figcnt=figcnt+l; 

figure  (  figcnt  ) 

loglog(r,dthetaed  ,  'o'  ,r,  dthetagd  ,  '  xk  '  ,r  ,dthetacn2ed  ,  '  *m'  ,r  ,dthetacn2gd  ,  '  .b  ') 
t  i  1 1  e  (  '  St  ructure  Function  of  Phase  comparison') 
x label (  ' r  (m)  '  ) 
ylabel  (  'Dtheta  '  ) 

legend  (  '  Edlen  '  ,  '  Gladstone— Dale  '  ,  'Cn~2  (Edlen )  '  ,  'Cn~2  (G— D)  '  ,  'Location  '  ,  'Bestoutside  ') 

%End  of  Routine 
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